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1  Research  Accomplishments 

SUMMARY: 

•  Developed  novel  multi-nesting,  time  splitting  and  implicit  relaxation  techniques  for  coupled 
mesoscale/microscale  ionospheric  dynamics;  resolved  secondary  Rayleigh- Taylor  instabilities. 

•  Characterized  neutral-plasma  interactions  and  turbulent  dynamics  of  plasma  bubbles/ionospheric 
layers  at  fine  scales. 

•  Created  methods  based  on  Lagrangian  coherent  structures  delineating  formation  of  scintillation 
producing  irregularities,  discovered  fine  scale  irregularity  patterns/plasma  pancakes  in  ionospheric 
flows;  described  non-Gaussian  spatially  varying  statistics  and  characterized  patchy  inhomogeneous 
dynamics  of  ionospheric  density  layers. 

•  Analyzed  stochastic  Lagrangian  dynamics  of  charged  flows  in  the  E-F  region  of  ionosphere; 
characterized  plasma  density  depletions  and  enhancements  induced  by  Lagrangian  dynamics,  and 
identified  barriers  for  ion/electron  transport. 

1.1  Multiscale  Modeling  and  Nested  Simulations  of  Three-Dimensional 

Ionospheric  Plasmas:  Rayleigh- Taylor  Turbulence  and  Non-Equilibrium 
Layer  Dynamics  at  Fine  Scales. 

Physics-based  predictive  modeling  and  novel  multi-nesting  computational  techniques  were  developed 
to  characterize  EM  wave  propagation  through  strongly  inhomogeneous  non-Kolmogorov  ionospheric 
media  (Mahalov  2014  [1]  and  references  therein).  Nested  numerical  simulations  of  ionospheric  plasma 
density  structures  associated  with  nonlinear  evolution  of  the  Rayleigh- Taylor  (RT)  instabilities  in 
Equatorial  Spread  F  (ESF)  were  carried  out.  The  high  resolution  in  targeted  regions  offered  by  the 
nested  model  was  able  to  resolve  scintillation  producing  ionospheric  irregularities  associated  with 
secondary  RT  instabilities  characterized  by  sharp  gradients  of  the  refractive  index  at  the  edges  of 
mixed  regions. 

For  the  limited  domain  and  nested  simulations,  the  lateral  boundary  conditions  were  treated  using 
novel  implicit  relaxation  techniques  applied  in  buffer  zones  where  the  density  of  charged  particles  for 
each  nest  was  relaxed  to  that  obtained  from  the  parent  domain.  Our  studies  also  focused  on  the 
charge- neutral  interactions  and  the  statistics  associated  with  stochastic  Lagrangian  dynamics.  We 
examined  the  organizing  mixing  patterns  for  plasma  flows  due  to  polarized  gravity  wave  excitations 
in  the  neutral  field,  using  Lagrangian  coherent  structures  (LCS).  LCS  objectively  depict  the  flow 
topology  and  the  extracted  scintillation  producing  irregularities  indicate  generation  of  ionospheric 
density  gradients,  due  to  accumulation  of  plasma. 

The  electron  density  in  the  ionosphere  varies  diurnally,  geographically,  seasonally,  with  sunspot 
number,  and  other  solar  phenomena.  The  total  electron  content  (TEC)  can  vary  by  two  orders  of 
magnitude  depending  on  the  time  and  location  of  observations.  Apart  from  the  variation  with  altitude, 
the  electron  density  varies  with  the  activity  level  of  the  sun,  time  of  the  year,  time  of  the  day,  and 
geographical  position.  In  addition,  a  number  of  stochastic  effects  and  scintillations  play  an  important 
role  in  the  influence  of  the  ionosphere  on  electromagnetic  wave  propagation. 

For  the  study  of  multi-scale  ionospheric  dynamics  over  a  limited  spatial  range  in  latitude,  longitude 
and  altitude,  nested  coupled  mesoscale/microscale  models  were  created  and  used  to  generate  high 
resolution  data  to  produce  accurate  predictions  and  forecasts.  We  developed  a  novel  multi-scale 
modeling  approach  or  nesting  to  capture  fine  scale  ionospheric  dynamics  where  models  that  are  able 
to  resolve  different  scales  are  coupled. 
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The  three-dimensional  equations  for  ionospheric  dynamics  included  coupling  of  neutral  fluid,  ion 
gas,  electron  gas,  and  electromagnetic  equations. 

3D  Navier-Stokes  Equations  for  Neutral  Fluid: 


^  +  v-u  =  o 

ot 

^+V.(Uti)  +  |^  +  (2nxU-JxB),  =  F„ 

^  +  V  ■  (Ut.)  +  +  (20  X  U  -  J  X  B),  =  Fv 

dW  dv 

^  +  V  •  (Urn)  +  ^  -  +  (2f7  X  U  -  J  X  =  Fhf 

dQ 

—  +  V  •  (U0)  =  Fe 
p  =  pKbT  /m  =  UKbT 
9  =  ^ 


In  these  equations,  U  =  {U,V,W)  is  the  coupled  velocity  (U  =  pu),  u  =  {u,v,w)  is  the  physical 
velocity  vector,  fl  is  the  earth  rotation,  J  is  the  density  current  of  charged  particles  (ions  and  electrons), 
B  is  the  earth  magnetic  field  vector,  T  is  the  temperature,  9  is  conserved  for  adiabatic  motion,  and 
0  is  its  coupled  variable  (0  =  p9),  p  is  the  pressure,  p  is  the  density,  g  is  the  acceleration  of  gravity 
vector,  m  is  the  mass,  n  is  the  number  density,  and  is  Boltzmanns  constant.  The  right-hand-side 
terms  F[/,  Fy,  Fw  and  Fg  represent  forcing  terms  arising  from  model  physics,  turbulent  mixing  and 
heating. 

The  above  equations  for  neutrals  are  coupled  with  equations  for  charged  fluid  which  include  ion 
and  electron  gases. 
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The  equations  for  the  ion  gas: 


dvi 


dni 

dt 


+  V  •  (njVj)  =  Pi 


Li 
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dTi 

dt 
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2  J_/,2^ 
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Qin  ^  ii  +  Q  ie 

Pi  TliKbTi. 


In  these  equations,  Vj  is  the  ion  velocity  vector,  E  is  the  electric  field  vector,  Tj  is  the  temperature, 
Pi  is  the  pressure,  n*  is  the  number  density  of  ion,  g  is  the  acceleration/gravity  vector,  Kj  is  the 
thermal  conductivity,  Vin  is  the  collision  frequency  of  ions  with  neutrals.  Here,  bg  =  ^  where  B  is  the 
geomagnetic  dipole  field  and  Bq  is  its  value  at  the  earths  surface,  and  s  is  the  coordinate  along  the 
magnetic  field  lines.  The  right-hand-side  terms  Pi  and  L,  represent  the  production  and  loss  coefficients 
of  ion  gas;  Qm,  Qn  and  Qie  represent  heating  terms  due  to  ion-neutral  collisions,  ion- ion  collisions, 
and  ion-electron  collisions  respectively. 

The  equations  for  the  electron  gas: 


Ue 


Pe-^  =  -nee{E  +  VeX  B)  +  UerUeg  -  Vpe  +  Heniel'enivi  “  Vg) 
dt  3  UeKb  ^  ds  \  ds 

Pe 

In  these  equations,  Ve  is  the  electron  velocity  vector,  E  is  the  electric  field  vector,  Tg  is  the 
temperature  of  electron  gas,  p^  is  the  pressure,  is  the  number  density  of  electron,  g  is  the  acceleration 
of  gravity  vector.  Kg  is  the  thermal  conductivity  for  electrons,  Uin  is  the  collision  frequency  of  electrons 
with  neutrals;  Qen  and  Q^i  represent  heating  terms  due  to  electron-neutral  collisions,  electron-ion 
collisions  respectively;  Qphe  is  due  to  photoelectron  heating. 

Electromagnetic  Equations: 


Qen  T  Q  ei  T  Qphe 


J  =  ^  UjqjYj  =  e  HiYi  -  UeVej 

V- J  =  0 

J  =  a  •  (E  u  X  B) 

E  =  -V0 

In  these  equations,  J  is  the  density  current,  E  is  the  electric  field  vector,  f  is  the  electric  potential, 
a  is  the  conductivity  tensor,  e  is  the  charge  of  the  electron,  and  qj  is  the  charge  of  the  j-th  species. 
The  electric  field  vector  (E)  and  the  potential  f  are  calculated  such  that  the  relation  V  •  J  =  0  is 
satisfied. 
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In  order  to  characterize  multi-scale  ionospheric  over  a  limited  spatial  range  in  latitude,  longitude 
and  altitude,  we  developed  nested  coupled  mesoscale/microscale  models  to  generate  high  resolution 
data  to  produce  accurate  predictions  and  forecasts.  Our  typical  limit  area  high  resolution  simulations 
focused  on  3D  ionospheric  domains  centered  on  geographical  regions  of  several  hundred  kilometers 
horizontally  and  encompassing  altitudes  from  80km  to  500km.  Nested  capabilities  allow  selected  time 
intervals  and  spatial  locations  to  run  at  higher  resolution.  The  smallest  mesoscale  domains  typically 
are  about  300km  300km  horizontally  in  the  ionosphere;  the  embedded  microscale  domains  are  100km 
100km  or  finer  and  zoomed  on  altitudes  encompassing  ionospheric  E-F  layers.  Running  global  models 
with  comparable  fine  resolution  to  resolve  ionospheric  meso/microscale  dynamics  in  targeted  regions 
is  prohibitive  due  to  high  computational  costs. 

Non-dimensional  numbers 

The  three-dimensional  physical  and  dynamical  processes  associated  with  polarized  neutral-plasma 
instabilities  were  analyzed.  The  horizontal  velocity  vector  (U{z),V{z))  (horizontally  and  time  av¬ 
eraged)  rotates  with  the  vertical  coordinate  (z).  Let  a{z)  denotes  the  angle  between  the  vector 

fhe  horizontal  wavevector.  The  following  expression  for  the  polarized  plasma 
Richardson  number  was  introduced: 


Here  the  definition  of  N  is  based  on  the  number  density  n  instead  of  temperature,  N'^  =  g  . 
The  polarized  plasma  Richardson  number  takes  into  account  horizontal  anisotropy  and  the  angle 
between  the  horizontal  wave  vectors  and  the  velocity  vectors  at  each  vertical  level.  In  the  case  of 
polarized  wind  fields  in  stably  stratified  neutral  environments,  the  polarized  Richardson  number  was 
rigorously  introduced  by  PI  in  earlier  publications.  The  above  expression  generalizes  this  definition 
to  ionospheric  flows,  [1],  where  electron  number  density  is  used  in  place  of  temperature  since  the 
ionosphere  is  stratified  with  respect  to  electron  density.  The  book  of  Kelley  [2],  ie  Page  159,  discusses 
a  plasma  Richardson  number  in  the  context  of  parallel  flows  in  the  ionosphere  and  the  role  of  velocity 
shear  in  triggering  convective  ionospheric  storms  (M.C.  Kelley,  Earths  Ionosphere,  Plasma  Physics 
and  Electrodynamics,  Amsterdam:  Elsevier,  2009).  The  above  expression  generalizes  his  definition. 

1.2  Neutral-Plasma  Interactions,  Ionospheric  Lagrangian  Coherent  Struc¬ 
tures  and  Turbulence  Statistics  of  Charged  Ionospheric  Layers. 

Three-dimensional  numerical  model  for  the  E-F  region  ionosphere  was  developed  and  used  to  study  the 
Lagrangian  dynamics  for  plasma  flows  in  this  region  (Mahalov  2014  [1]  and  references  therein).  These 
studies  focused  on  charge-neutral  interactions  and  the  statistics  associated  with  stochastic  Lagrangian 
motion.  In  particular,  the  organizing  mixing  patterns  for  plasma  flows  due  to  polarized  gravity 
wave  excitations  in  the  neutral  field  were  delineated,  using  Lagrangian  coherent  structures  (LCS) 
techniques.  LCS  objectively  depict  the  flow  topology  and  the  extracted  attractors  indicate  generation 
of  ionospheric  density  gradients,  due  to  accumulation  of  plasma.  Using  Lagrangian  measures  such  as 
the  finite-time  Lyapunov  exponents,  the  Lagrangian  skeletons  for  mixing  in  plasma  were  delineated, 
and  charged  fronts  were  characterized  (strong  scintillation  regions  in  space).  With  polarized  neutral 
wind,  we  found  that  the  corresponding  plasma  velocity  is  also  polarized.  Moreover,  the  polarized 
neutral  winds,  coupled  with  stochastic  Lagrangian  motion,  give  rise  to  polarized  density  fronts  in 
plasma.  Statistics  of  these  trajectories  showed  high  level  of  non-Gaussianity.  This  includes  clear 
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signatures  of  variance,  skewness,  and  kurtosis  of  displacements  taking  polarized  structures  aligned 
with  the  polarized  ionospheric  density  fronts,  the  latter  being  strongly  anisotropic. 

Three-dimensional  numerical  studies  of  the  E  and  lower  F  region  ionosphere  coupled  with  the 
neutral  atmosphere  dynamics  were  conducted.  A  model  was  developed  to  resolve  the  transport  pat¬ 
terns  of  plasma  density  coupled  with  neutral  atmospheric  dynamics.  Inclusion  of  neutral  dynamics 
in  the  model  allowed  to  examine  the  charge-neutral  interactions  over  the  full  non-equilibrium  evo¬ 
lution  cycle  of  an  inertial  gravity  wave  when  the  background  flow  spins  up  from  rest,  saturates  and 
eventually  breaks.  Using  Lagrangian  analyses,  the  mixing  patterns  of  the  ionospheric  responses  and 
the  formation  of  ionospheric  layers  were  delineated.  The  corresponding  plasma  density  in  this  flow 
develops  complex  wave  structures  and  fine  scale  patches  of  scintillation  producing  irregularities  during 
the  gravity  wave  breaking  event. 

Additional  details  can  be  found  in  the  Appendix,  Mahalov  2014  [1]  and  references  therein. 
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21,  042901. 
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and  zonal  flows,  Europ.  Journal  of  Mechanics  B/Fluids,  vol.  37,  Pages  48-58. 

•  Durazo,  J.,  E.  Kostelich,  A.  Mahalov  and  W.Tang.  2015.  Assessing  a  local  ensemble  trans¬ 
form  Kalman  Alter:  Observing  system  experiments  with  an  ionospheric  model,  special  volume 
on  Geophysical  and  Astrophysical  turbulence,  accepted  for  publication,  Cambridge  University 
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3  Selected  Invited  and  Plenary  Presentations  at  Scientific 
Meetings  and  Conferences 

•  12th  Conference  on  Space  Weather,  95th  Annual  Meeting  of  the  American  Meteorological  Soci¬ 
ety,  Phoenix  AZ,  4-8  January,  2015  (to  be  presented) 

•  Institute  for  Pure  and  Applied  Mathematics  (IPAM),  UCLA,  Geophysical  and  Astrophysical 
Turbulence  Workshop,  several  invited  presentations  in  the  IPAM  Fall  2014  Turbulence  Program 

•  International  Center  for  Theoretical  Physics,  Mixing  in  Rapidly  Changing  Environments:  Prob¬ 
ing  Matter  at  the  Extremes,  Trieste,  Italy,  4-9  August,  2014 

•  Optical  Society  of  America  Congress  on  Imaging  and  Applied  Optics,  Propagation  through  and 
Characterization  of  Distributed  Volume  Turbulence  Symposium,  Seattle,  13-17  July,  2014 

•  Advances  in  Mathematical  Eluid  Mechanics:  Stochastic  and  Deterministic  Methods,  Portugal 
2014 

•  AFOSR  Space  Science  Program  Annual  Meeting,  13-14  January,  2014 

•  Annual  Meeting  of  the  Division  of  Plasma  Physics  of  the  American  Physical  Society,  Denver, 
November  11-15,  2013 

•  Optical  Society  of  America  Congress  on  Imaging  and  Applied  Optics,  Arlington,  23-27  June, 
2013 

•  Center  for  Sci.  Comp,  and  Math.  Modeling  Colloquium,  Univ.  of  Maryland,  March  13,  2013 

•  DoD  National  Reconnaissance  Office  (NRO)  Technology  Seminar  Series,  Washington  DC,  March 
21,  2013 

•  Scientific  Computing  Colloquium,  Brown  University,  November  19,  2012 

4  Personnel 

•  Erwin  Suazo,  Research  Scientist.  Currently  tenure-track  Assistant  Professor  at  the  University 
of  Puerto  Rico,  Mayaguez.  Erwin  established  scientific  contacts  with  the  Arecibo  observatory 
group  at  the  National  Astronomy  and  Ionosphere  Center,  Puerto  Rico 

•  Juan  Durazo,  PhD  student,  MS  completed  in  2013,  PhD  thesis  defense  expected  in  Spring  2015. 
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Abstract 

Multiscale  modeling  and  high  resolution  three-dimensional  simulations  of  nonequihbrium 
ionospheric  dynamics  are  major  frontiers  in  the  field  of  space  sciences.  The  latest  developments  in 
fast  computational  algorithms  and  novel  numerical  methods  have  advanced  reliable  forecasting  of 
ionospheric  environments  at  fine  scales.  These  new  capabilities  include  improved  physics-based 
predictive  modeling,  nesting  and  imphcit  relaxation  techniques  that  are  designed  to  integrate 
models  of  disparate  scales.  A  range  of  scales,  from  mesoscale  to  ionospheric  microscale,  are 
included  in  a  3D  modeling  framework.  Analyses  and  simulations  of  primary  and  secondary 
Rayleigh-Taylor  instabilities  in  the  equatorial  spread  F  (ESF),  the  response  of  the  plasma  density 
to  the  neutral  turbulent  dynamics,  and  wave  breaking  in  the  lower  region  of  the  ionosphere  and 
nonequilibrium  layer  dynamics  at  fine  scales  are  presented  for  coupled  systems  (ions,  electrons  and 
neutral  winds),  thus  enabling  studies  of  mesoscale/micro scale  dynamics  for  a  range  of  altitudes 
that  encompass  the  ionospheric  E  and  F  layers.  We  examine  the  organizing  mixing  patterns  for 
plasma  flows,  which  occur  due  to  polarized  gravity  wave  excitations  in  the  neutral  field,  using 
Lagrangian  coherent  structures  (LCS).  LCS  objectively  depict  the  flow  topology  and  the  extracted 
scintillation-producing  irregularities  that  indicate  a  generation  of  ionospheric  density  gradients,  due 
to  the  accumulation  of  plasma.  The  scintillation  effects  in  propagation,  through  strongly 
inhomogeneous  ionospheric  media,  are  induced  by  trapping  electromagnetic  (EM)  waves  in 
parabolic  cavities,  which  are  created  by  the  refractive  index  gradients  along  the  propagation  paths. 

Keywords:  plasma  dynamics,  ionospheric  layers,  Rayleigh-Taylor  instabilities,  nested 
simulations 


1.  Introduction 

Turbulent  hydrodynamic  mixing,  induced  by  Rayleigh- 
Taylor  (RT)  instabilities,  occurs  in  settings  that  are  as  varied 
as  exploding  stars  (supemovae),  inertial  confinement  fusion 
(ICE)  and  macroscopic  flows  in  fluid  dynamics,  such  as 
ionospheric  plasmas  ([1^,  11,  58,  95,  99,  100,  122]).  Since 

^  The  Wilhoit  Foundation  Dean’s  Distinguished  Professor. 


the  discovery  of  the  plasma  instability  phenomenon  that 
occurs  in  the  nighttime  equatorial  F-region  ionosphere,  which 
is  revealed  by  rising  plumes  that  are  identified  as  large-scale 
depletions  or  bubbles,  considerable  efforts  have  been  made  in 
the  development  of  computer  models  that  simulate  the  gen¬ 
eration  and  evolution  of  the  Equatorial  Spread  F  (ESF) 
dynamics  ([9,  11,  17,  29,  30,  41,  49,  58,  75,  76,  83-85,  87, 
89,  99,  100,  104,  110]).  The  book  by  M  C  Kelley  [50], 
contains  a  comprehensive  list  of  references  to  2009. 
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The  ionosphere  is  a  dynamic  mixture  of  ions,  electrons 
and  neutral  gases,  surrounding  the  Earth  in  the  altitude  range 
of  approximately  90  km  to  beyond  1000  km.  The  ionosphere 
can  also  be  viewed  as  a  transition  region  from  the  Earth’s 
lower  atmospheric  regions  (i.e.,  the  troposphere,  stratosphere 
and  mesosphere)  to  the  outer  space  environment  (i.e.,  the 
magnetosphere).  As  such,  the  ionosphere  acts  to  mediate  and 
transmit  external  forces  and  drivers  from  below  and  above. 
Understanding  and  modeling  the  ionosphere  is  a  daunting  but 
important  challenge  because  of  its  critical  role  in  space  sci¬ 
ences.  Significant  progress  has  been  made  in  the  development 
of  computational  models  of  the  ionosphere  over  the  past 
several  decades.  However,  the  current  modeling  capabilities 
need  further  development,  especially  with  regard  to  coupling 
the  physical  processes  that  occur  over  a  large  range  of  tem¬ 
poral  and  spatial  scales  that  characterize  the  ionosphere.  In 
particular,  advances  are  needed  in  the  ability  to  make  accurate 
high  resolution  forecasts  over  limited-area  ionospheric 
environments,  which  are  required  for  operational  commu¬ 
nication,  navigation  and  imaging  applications. 

The  ionosphere  involves  interactions  between  phenom¬ 
ena  of  varying  scale  sizes  ([6,  12,  19-21,  26,  33-34,  36-41, 
46,  48,  52,  53,  72,  77-90,  96-122]).  Large-scale  variations, 
like  the  solar  cycle,  seasonal  effects  and  tidal  effects,  all  drive 
large-scale  changes  in  the  global  structure.  Such  changes 
define  the  mean  state  (climate)  of  the  ionospheric  system. 
These  processes  have  characteristic  spatial  scales  greater  than 
a  thousand  kilometers  and  time  scales  that  range  from  several 
hours  to  a  few  years.  General  circulation  models  have  been 
developed  to  simulate  these  large-scale  structures  (e.g., 
[14,  18,  22,  23,  47,  48  86]).  These  global  models,  together 
with  the  large  observational  data  sets  that  have  been  accu¬ 
mulated  over  the  years,  have  led  to  a  much  greater  under¬ 
standing  of  large-scale  structures  in  the  ionosphere  and  to  the 
response  of  these  structures  to  variations  in  geophysical 
inputs.  Space  weather  is  the  perturbation  of  the  ionosphere 
and  the  thermosphere  from  its  long-term  global  mean  state. 
These  perturbations  involve  not  only  large-scale  variations, 
but  also  mesoscale  and  small-scale  processes  that  occur 
locally  and  may  have  short  periods.  Mesoscale  and  small- 
scale  processes,  such  as  auroral  arcs  and  the  mid-latitude 
electron  density  trough,  affect  not  only  local  plasma  and 
neutral  distributions  but  also  large-scale  structures  through 
dynamical  and  energetic  coupling.  Such  couplings  between 
processes  of  different  scales  occur  in  the  equatorial  region  and 
in  the  mid-latitudes  and  high  latitudes.  The  dynamics  of  large 
and  mesoscale  neutral  winds  are  often  characterized  by  the 
presence  of  zonal  jets  and  anisotropic  turbulence,  i.e.  [24], 
and  [13].  Improving  our  understanding  of  neutral/plasma 
interactions  at  ionospheric  meso/micro scales  is  a  challenge 
for  ionospheric  science  studies.  M  C  Kelley  ([50],  p.  456) 
states,  ‘we  simply  do  not  understand  neutral/plasma  interac¬ 
tions  on  our  own  planet.’  Current  models  are  limited  because 
of  their  use  of  empirical  models  of  the  ionosphere  for  pre¬ 
scribed  neutral  conditions,  which  decouple  the  dynamics  of 
neutral  and  charged  particles.  They  are  also  limited  by  the 
assumption  of  neglecting  the  inertial  terms  and  by  imposing 
isothermal  conditions  for  both  charged  and  neutral  species. 


The  electron  density  in  the  ionosphere  varies  diumally, 
geographically,  seasonally  and  with  the  Sunspot  number, 
along  with  other  solar  phenomena.  The  total  electron  content 
(TEC)  can  vary  by  two  orders  of  magnitude,  depending  on  the 
time  and  the  location  of  the  observations.  Apart  from  the 
variation  with  altitude,  the  electron  density  varies  with  the 
activity  level  of  the  Sun,  the  time  of  the  year,  the  time  of  the 
day  and  the  geographical  position.  In  addition,  a  number  of 
stochastic  effects  and  scintillations  play  an  important  role  in 
the  influence  of  the  ionosphere  on  electromagnetic  wave 
propagation.  At  low  latitude,  an  important  scintillation  source 
is  E-spread.  E-spread  is  caused  by  rod-shaped  magnetic  field- 
aligned  bubbles,  which  are  formed  in  the  E-layer  just  after 
sunset  and  have  a  lifetime  of  2-3  h.  The  edges  of  the  E-spread 
bubbles  are  highly  unstable  and  can  be  the  source  of  intensity 
scintillations.  E-spread  is  more  prevalent  during  equinoxes 
and  summers,  occurs  preferentially  during  magnetically  quiet 
periods  and  increases  with  increasing  sun  activity.  At  high 
latitude,  the  aurora  can  cause  severe  distortions  of  the  elec¬ 
tromagnetic  (EM)  waveforms.  The  aurora  is  the  result  of 
high-energy  electrons  from  the  solar  wind,  which,  at  polar 
regions,  can  sometimes  break  through  the  barrier  of  the 
Earth’s  magnetic  field.  These  electrons  ionize  atoms  and 
therefore  cause  the  electron  density  to  increase. 

Eor  the  study  of  multiscale  ionospheric  dynamics  over  a 
limited  spatial  range  in  latitude,  longitude  and  altitude,  nested 
coupled  mesoscale/microscale  models  need  to  be  created  and 
used  in  order  to  generate  high  resolution  data,  which  is  then 
used  to  produce  accurate  predictions  and  forecasts.  The 
challenge  of  multiscale  nesting  and  High  Performance  Com¬ 
puting  (HPC)  simulations  for  the  dynamically  evolving  lim¬ 
ited-area  ionospheric  environments  is  to  achieve  the  robust 
predictability  and  ensemble  forecasting  of  high-impact  iono¬ 
spheric  events.  Eurthermore,  these  high  resolution  regions 
need  to  be  resolved  even  further  by  coupling  to  advanced 
fluid,  hybrid,  and/or  Lagrangian  particle  tracking  codes, 
which  provide  the  necessary  physics  to  study  ionospheric 
processes  at  fine  scales.  This  is  a  daunting  but  important 
challenge  to  meet.  Significant  advances  in  the  computation  of 
atmospheric  and  environmental  flows  at  lower  atmospheric 
levels  have  been  achieved  during  the  last  decade.  A  dramatic 
increase  in  computer  power  has  facilitated  the  development  of 
numerical  codes  that  have  the  capability  to  resolve  small-scale 
atmospheric  processes.  This  was  achieved  by  the  imple¬ 
mentation  of  fast  computational  algorithms  and  nesting 
techniques  with  multiple  domains  that  resolve  fine  horizontal 
and  vertical  scales  and  by  the  improvement  of  physics-based 
predictive  modeling  for  shear- stratified  flows  in  the  upper 
troposphere  and  lower  stratosphere  ([43,  59,  60]).  Multiscale 
modeling  of  the  Earth’s  ionosphere  has,  however,  received 
less  attention  compared  to  lower  atmospheric  levels,  and  it 
remains  one  of  the  major  frontiers  of  space  sciences.  One  of 
the  main  difficulties  researchers  face  in  the  development  of  a 
nested  modeling  approach  for  ionospheric  equations  is  the 
specification  of  the  boundary  conditions.  Additional  compli¬ 
cations  in  the  computations  of  ionospheric  flows  arise  from 
the  equation  for  the  electric  potential,  which  is  a  diagnostic 
variable  computed  from  nonlinear  elliptic  equations. 
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In  this  review,  we  describe  a  novel  multiscale  modeling 
approach,  or  nesting,  to  capture  fine-scale  ionospheric 
dynamics,  in  which  models  that  are  able  to  resolve  different 
scales  are  coupled.  The  prognostic  fields  at  the  boundaries  of 
the  nested  model  are  specified  from  a  large-scale  model.  The 
large  domain  model,  with  a  coarse  resolution,  is  used  to 
predict  large-scale  structures,  while  a  limited  domain  (or 
nested)  model  with  boundary  conditions  interpolated  from  the 
coarse  grid  (or  parent)  model  is  used  over  smaller  domains 
with  finer  resolutions.  Small-scale  structures,  which  are  not 
resolved  in  the  coarse  grid  model  and  therefore  need  to  be 
represented  by  some  kind  of  parameterization,  are  explicitly 
resolved  in  the  nested  model;  this  is  the  improvement 
achieved  by  this  approach.  This  methodology  was  applied, 
with  success,  to  the  multiscale  resolution  of  the  atmospheric 
flows  in  the  upper  troposphere  and  lower  stratosphere 
([43,  59-66]).  The  extension  of  multi-nested  computational 
methods  to  the  ionosphere  presents  significant  challenges,  as 
the  presence  of  ions  and  electrons  increases  the  complexity  of 
the  model  equations. 

Here,  we  present  high  resolution  multiscale  nested 
simulations  of  ionospheric  dynamics,  in  which  nested  fields 
are  relaxed  toward  large  domain  fields.  We  describe  a  new 
multiscale  computational  method  based  on  the  flow  relaxation 
scheme,  in  which  the  relaxation  is  implemented  as  an  implicit 
correction  ([59]).  We  show  examples  that  demonstrate  that 
this  method  is  very  effective  and  robust  in  the  context  of 
ionospheric  flows.  Computational  results  in  real  ionospheric 
Equatorial  Spread  F  (ESF)  conditions  are  presented,  which 
demonstrate  the  ability  of  the  proposed  nested  approach  to 
resolve  the  multiscale  physics  of  strongly  nonlinear  structures 
observed  in  the  ESF  and  to  improve  the  resolution  of  the 
primary  and  secondary  ionospheric  Rayleigh-Taylor  (RT) 
instabilities  and  non-equilibrium  plasma  dynamics 
([58,  99,  100]).  We  examine  the  organizing  mixing  patterns 
for  the  plasma  flows  due  to  polarized  gravity  wave  excitations 
in  the  neutral  field,  using  Lagrangian  coherent  structures 
(ECS).  ECS  objectively  depict  the  flow  topology,  and  the 
extracted  scintillation  producing  irregularities  indicate  a 
generation  of  ionospheric  density  gradients  due  to  the  accu¬ 
mulation  of  plasma.  The  scintillation  effects  are  induced  by 
trapping  electromagnetic  (EM)  waves  in  parabolic  cavities, 
created  by  the  refractive  index  gradients  along  the  propaga¬ 
tion  paths  ([54,  57,  68]). 

This  paper  is  organized  as  follows.  The  model  formula¬ 
tion  and  the  computational  approach  are  presented  in 
section  2.  The  performance  and  the  ability  of  the  nested 
simulations  to  resolve  small-scale  structures  associated  with 
Rayleigh-Taylor  instabilities  and  the  ESF  are  demonstrated  in 
section  3.  In  section  4,  our  studies  focus  on  the  charge-neutral 
interactions  and  the  associated  stochastic  Lagrangian 
dynamics.  In  particular,  we  examine  the  organizing  mixing 
patterns  for  plasma  flows  due  to  polarized  gravity  wave 
excitations  in  the  neutral  field,  using  Lagrangian  coherent 
structures  (ECS)  methods.  Finally,  the  summary  and  chal¬ 
lenges  are  discussed  in  section  5. 


2.  Multiscale  coupled  mesoscale/microscale  limited- 
area  ionospheric  models:  nesting  and  implicit 
relaxation  computational  techniques 


2.1.  Governing  equations 

The  three-dimensional  equations  for  ionospheric  dynamics 
include  the  coupling  of  the  neutral  fluid,  the  ion  gas,  the 
electron  gas  and  the  electromagnetic  equations  ([50,  88]). 

2. 1. 1.  Equations  for  the  neutrai  fiuid. 

^  +  V  .  u  =  0 
dt 

—  +V-  {Uu)  +  -^  +  (2X3  X  t7  -  /  X  b)  =  Fy 

dt  dx  X 

—  +  V  ■{ijv)  +  —  +  {2QxU  -J  xb)  =  Fy 

dt  ^  ^  dy  ^  ’y 

—  +V-  (Uw)  +  ^  -  pg  +  {2QxU  -J  xb]  =F^ 

dt  ^  ^  dz  ^ 

f +F.(™)  =  F, 

p  =  pk^T/m  =  nkf^T 

e  =  ^. 

p' 

In  these  equations,  U  =  (U,V,W)  is  the  coupled  velocity 
(U  =  pu),u  =  (u,v,w)  is  the  physical  velocity  vector,  Q  is  the 
Earth’s  rotation,  J  is  the  density  current  of  the  charged 
particles  (ions  and  electrons),  B  is  the  Earth  magnetic  field 
vector,  T  is  the  temperature,  6  is  conserved  for  adiabatic 
motion  and  0  is  its  coupled  variable  {0  =  p6).  p  is  the 
pressure,  p  is  the  density,  g  is  the  acceleration  of  gravity 
vector,  m  is  the  mass,  n  is  the  number  density  and  is 
Boltzmann’s  constant.  The  right-hand  side  terms  Fy, 
and  F^  represent  forcing  terms  that  arise  from  model  physics, 
turbulent  mixing  and  heating. 

The  above  equations  for  neutrals  are  coupled  with 
equations  for  the  charged  fiuid,  which  include  ion  and 
electron  gases. 


2. 1.2.  Equations  for  the  ion  gas. 


dn.  ^ 

-^+V  ■  (n^v,)  =  P>-L, 

Pi~~  —  (f  +  X  b)  +  n-nijg  —  Vp^  +  [u  —  it) 


dX  ^  ^  2  - 

—  +  it  •  VTF  -TV 
dt  '  3  ' 

=  Qin  +  Qii  +  Qie 

Pi  = 


In  these  equations,  if  is  the  ion  velocity  vector,  E  is  the 
electric  field  vector,  Ti  is  the  temperature,  pt  is  the  pressure, 
n.  is  the  number  density  of  ion,  g  is  the  acceleration  of  the 
gravity  vector,  /c.  is  the  thermal  conductivity  and  is  the 
collision  frequency  of  the  ions  with  the  neutrals.  Here, 
where  B  is  the  geomagnetic  dipole  field,  and  is 
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Figure  1 .  Horizontal  and  vertical  cells  of  the  C  grid  staggering. 


its  value  at  the  Earth’s  surface,  and  s  is  the  coordinate  along 
the  magnetic  field  lines.  The  right-hand  side  terms  and  L. 
represent  the  production  and  loss  coefficients  of  the  ion  gas. 

Q..  and  represent  the  heating  terms,  which  are  due  to 
ion-neutral  collisions,  ion-ion  collisions  and  ion-electron 
collisions,  respectively. 


2.1.3.  Equations  for  the  electron  gas. 


=  Z"- 


dv  ^  ^  ^  ^ 

=  -n,e{E  +  V,  X  B)  +  n^m^g 

-  Vp^  +  -  i?) 


phe 


B 


In  these  equations,  is  the  electron  velocity  vector,  E  is  the 
electric  field  vector,  is  the  temperature  of  the  electron  gas, 
Pe  is  the  pressure,  is  the  number  density  of  the  electron,  g  is 
the  acceleration  of  the  gravity  vector,  k;  is  the  thermal 
conductivity  for  electrons  and  is  the  collision  frequency  of 
the  electrons  with  the  neutrals.  and  2^.  represent  the 
heating  terms,  which  are  due  to  the  electron-neutral  collisions 
and  to  the  electron-ion  collisions,  respectively.  is  due  to 
photoelectron  heating. 


2.1.4.  Electromagnetic  equations. 


^  ^ 

j 


V  i 


F  •  7  =  0 

J  =  (79  (^E-\-uxB^ 
E  =  —  V(p. 


In  these  equations,  7  is  the  density  current,  E  is  the  electric 
field  vector,  cf)  is  the  electric  potential,  o  is  the  conductivity 
tensor,  e  is  the  charge  of  the  electron  and  qj  is  the  charge  of 


the  7th  species.  The  electric  field  vector  (E)  and  the  potential 
(f)  are  calculated  so  that  the  relation  V  •  7  =  0  is  satisfied. 

In  order  to  characterize  multiscale  ionospheric  dynamics 
over  a  limited  spatial  range  in  the  latitude,  longitude  and 
altitude,  we  use  nested  coupled  mesoscale/microscale  models 
to  generate  high  resolution  data,  which  is  used  to  produce 
accurate  predictions  and  forecasts.  Our  typical  limit  area  high 
resolution  simulations  focus  on  a  3D  ionospheric  domain, 
which  is  centered  on  a  geographical  region  that  spans  several 
hundred  kilometers  horizontally  and  encompasses  altitudes 
from  80  km  to  500  km.  Nested  capabilities  allow  selected 
time  intervals  and  spatial  locations  to  run  at  a  higher 
resolution.  The  smallest  mesoscale  domains  are  typically 
about  300  km  X  300  km  horizontally  in  the  ionosphere;  the 
embedded  microscale  domains  are  100  km  x  100  km  or  finer 
and  are  zoomed  on  altitudes  that  encompass  the  ionospheric 
E-F  layers.  Running  global  models  with  comparable  fine 
resolution  in  order  to  resolve  ionospheric  meso/micro scale 
dynamics  in  targeted  regions  is  prohibitive  due  to  the 
computational  costs. 


2.2.  Model  numerics 

The  solver  for  neutral  and  charged  gases  uses  a  time-split 
integration  scheme.  Slow  or  low-frequency  modes  are  inte¬ 
grated  using  a  third-order  Runge  Kutta  (RK3)  time  integration 
scheme,  while  the  high-frequency  modes  are  integrated  over 
smaller  time  steps  in  order  to  maintain  numerical  stability. 
The  horizontally  propagating  acoustic  modes  and  fast  waves 
are  integrated  using  a  forward-backward  time  integration 
scheme,  and  vertically  propagating  acoustic  modes  are  inte¬ 
grated  using  a  vertically  implicit  scheme  (using  the  acoustic 
time  step).  This  time- splitting  is  similar  to  what  is  described 
in  the  atmospheric  codes  and  analyzed  in  [60,  92].  The 
acoustic-mode  integration  is  cast  in  the  form  of  a  correction  to 
the  RK3  integration.  The  RK3  scheme  integrates  a  set  of 
partial  differential  equations  by  using  a  predictor-corrector 
formulation.  The  high-frequency  acoustic  modes  would 
severely  limit  the  RK3  time  step.  We  use  the  time-split 
approach  in  order  to  circumvent  this  time  step  limitation. 
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Additionally,  to  increase  the  accuracy  of  the  splitting,  we 
integrate  a  perturbation  form  of  the  governing  equations  by 
using  smaller  acoustic  time  steps  within  the  RK3  large-time- 
step  sequence.  To  form  the  perturbation  equations  for  the 
RK3  time-split  acoustic  integration,  we  define  small  time  step 
variables  that  are  deviations  from  the  most  recent  RK3  pre¬ 
dictor  {y/"  =  y/  —  y/^). 

The  spatial  discretization  in  the  solver  uses  C  grid  stag¬ 
gering  for  the  variables,  as  shown  in  figure  1.  That  is,  normal 
velocities  are  staggered  one-half  grid  length  from  the  ther¬ 
modynamic  variables.  The  variable  indices  (/,  7,  k)  indicate 
variable  locations.  We  denote  the  points  where  the  density  is 
located  as  thermodynamic  points,  and  likewise,  we  denote  the 
locations  where  u,  v,  and  w  are  defined  as  u  points,  v  points 
and  w  points,  respectively.  The  mass  and  the  number  density 
for  neutrals,  ions  and  electrons  are  defined  at  the  (/,  7,  k) 
points.  The  pressure  p  is  computed  at  the  thermodynamic 
points.  The  grid  spacing  Ax  and  Ay  are  constants  in  the  model 
formulation.  The  vertical  grid  spacing  Az  is  not  a  fixed  con¬ 
stant;  it  is  specified  in  the  initialization. 

The  RK3  advection  uses  a  third  order  in  the  vertical 
accurate  spatial  discretizations  and  a  fifth  order  in  the  hor¬ 
izontal  accurate  spatial  discretizations  of  the  flux  divergence 
for  momentum,  scalars  and  temperature,  using  the  RK3  time- 
integration  scheme.  The  formulation  of  the  odd-order 
schemes  are  comprised  of  the  next  higher  (even)  order  cen¬ 
tered  scheme,  plus  an  upwind  term  that,  for  a  constant 
transport  mass  flux,  is  a  diffusion  term  of  that  next  higher 
(even)  order  with  a  hyper- viscosity  that  is  proportional  to  the 
Courant  number.  The  diffusion  term  is  the  leading  order  error 
term  in  the  flux  divergence  discretization.  The  Positive- 
Definite  Limiter  for  the  RK3  advection  of  charged  and  neutral 
chemical  species  or  other  tracer  species  should  remain  posi¬ 
tive  definite;  that  is,  negative  masses  should  not  be  permitted. 
The  Runge-Kutta  transport  integration,  defined  by  the  time 
stepping  algorithm  and  combined  with  the  flux  divergence 
operator,  is  conservative,  but  it  does  not  guarantee  positive 
definiteness.  Any  negative  values  are  offset  by  the  positive 
mass,  such  that  the  mass  is  conserved.  In  many  physics 
options,  negative  concentrations  and  densities  are  set  to  zero, 
which  results  in  an  increase  in  the  mass  of  that  species.  A 
positive-definite  flux  renormalization,  applied  on  the  final 
Runge-Kutta  transport  step,  can  be  used  to  remove  this 
unphysical  effect  from  the  RK3  scalar  transport  scheme.  The 
model  is  typically  integrated  with  a  fixed  time  step,  which  is 
chosen  to  produce  a  stable  integration.  During  any  time  in  the 
integration,  the  maximum  stable  time  step  is  likely  to  be 
larger  than  the  fixed  time  step.  An  adaptive  time  stepping 
capability  is  introduced  so  that  the  RK3  time  step  can  be 
chosen  based  on  the  temporally  evolving  wind  fields.  The 
adaptively  chosen  time  step  is  usually  larger  than  the  typical 
fixed  time  step;  hence,  the  dynamics  integrate  faster,  and  the 
physics  parameterizations  are  called  less  often.  Also,  the  time- 
to-completion  of  the  simulation  can  be  substantially  reduced. 

Recently,  a  fast  time  stepping  method,  based  on  leapfrog 
and  a  new  high  order  implicit  time  filter,  was  developed  in 
[73].  The  amplitude  errors  in  this  scheme  are  of  the  third 
order,  which  is  comparable  to  the  amplitude  errors  of  the  third 


order  RK  scheme.  The  new  scheme,  however,  only  requires 
one  function  evaluation  per  time  step  as  opposed  to  the  RK3 
time  step,  which  is  more  expensive  since  it  requires  three 
evaluations  for  each  time  step  (the  new  scheme  is  three  times 
faster  than  the  currently  used  schemes  and  incorporates  the 
same  level  of  accuracy).  This  makes  the  numerical  method  in 
[73]  a  potential  alternative  for  applications  in  high  resolution 
computational  modeling  and  in  the  forecasting  of  ionospheric 
dynamics. 

2.3.  Multi-Nesting 

In  the  study  of  limited-area  ionospheric  environments,  one 
key  challenge  is  the  development  of  nesting  methods  that  use 
robust  computational  techniques  to  control  numerical  errors, 
such  as  the  ones  that  are  generated  at  the  boundaries  of  the 
nested  models.  These  errors  are  inevitable  in  any  nested  or 
downscaled  simulation  because  the  coarse  grid  fields  that  are 
specified  at  the  boundaries  of  the  nested  grid  are  not  con¬ 
sistent  with  the  finer  fields  that  are  simulated  by  the  nest.  If 
not  treated  effectively,  these  errors  will  propagate  into  the 
interior  of  the  nested  domain,  thereby  resulting  in  the  poor 
quality  of  the  numerical  simulations.  These  errors  are  parti¬ 
cularly  important  when  instabilities  dominate  the  dynamics 
near  the  boundaries  of  the  nested  domain.  Thus,  numerical 
methods  that  control  and  reduce  the  propagation  of  these 
errors  are  essential  for  reliable  high  resolution  regional  pre¬ 
dictions  in  the  ionosphere. 

There  are  two  main  techniques  that  are  used  in  atmo¬ 
spheric  and  oceanic  downscaling  models  to  improve  resolu¬ 
tion  over  limited  areas.  In  dynamically  adaptive  methods,  the 
spatial  resolution  constantly  changes  with  time  by  coarsening 
or  refining  the  grid  spacing,  depending  on  the  local  conditions 
([8,  15]).  The  adaptive  methods  are  not  well-established  in  the 
atmospheric  modeling  systems  for  several  reasons  ([8]):  (i) 
adaptive  techniques  can  incur  massive  overhead  due  to 
indirect  data  addressing  and  to  the  additional  efforts  needed 
for  grid  handling,  which  increase  the  cost  of  the  real-time 
forecasting  or  the  long-term  predictions;  (ii)  physical  para¬ 
meterizations  of  subgrid  processes  are  usually  optimized  for  a 
specific  grid  resolution,  which  makes  it  difficult  to  use 
dynamically  and  temporally  refined  or  coarsened  grids.  An 
important  class  of  numerical  methods  uses  nesting  to  improve 
spatial  resolution  over  a  limited  area.  Nesting  techniques  are 
widely  used  in  atmospheric  ([42,  59,  60,  64,  66,  92-94])  and 
oceanic  ([28,  91])  models.  Large  domain  models  with  coarse 
resolution  are  used  to  predict  large-scale  dynamics,  while 
limited-area  models  with  boundary  conditions  that  are  inter¬ 
polated  from  coarse  grids  are  used  over  embedded  domains 
with  a  finer  resolution.  Small-scale  processes,  which  are  not 
resolved  in  a  coarse  grid  model  and  therefore  need  to  be 
represented  by  using  subgrid-scale  parameterizations,  may  be 
explicitly  resolved  in  the  nested  model;  this  is  the  improve¬ 
ment  achieved  by  nesting  techniques. 

The  horizontal  and  vertical  nesting  capabilities  allow  the 
resolution  to  be  focused  over  a  region  of  interest  in  the 
ionosphere  by  introducing  an  additional  grid  (or  grids)  into 
the  simulation.  Vertical  and  horizontal  refinements  are 
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Figure  2.  (a)  Horizontal  Arakawa-C  grid  staggering  for  a  portion  of  a  parent  domain  and  an  imbedded  nest  domain  with  a  3:1  grid  size  ratio. 
The  solid  lines  denote  coarse  grid  cell  boundaries,  and  the  dashed  lines  are  the  boundaries  for  each  fine  grid  cell.  The  horizontal  components 
of  velocity  {U  and  V)  are  defined  along  the  normal  cell  face,  and  the  thermodynamic  variables  n  are  defined  at  the  center  of  the  grid  cell  (each 
square).  The  variables  along  the  interface  between  the  coarse  and  the  fine  grid  define  the  locations  where  the  specified  lateral  boundaries  for 
the  nest  are  in  effect,  (b)  Vertical  Arakawa-C  grid  staggering  for  a  portion  of  a  parent  domain  and  an  imbedded  nest  domain  with  a  3:1  grid 
size  ratio.  The  solid  lines  denote  coarse  grid  cell  boundaries,  and  the  dashed  lines  are  the  boundaries  for  each  fine  grid  cell.  The  eastward  and 
the  vertical  components  of  velocity  {U  and  W)  are  defined  along  the  normal  cell  face,  and  the  thermodynamic  variables  n  are  defined  at  the 
center  of  the  grid  cell  (each  square).  The  variables  along  the  interface  between  the  coarse  and  the  fine  grid  define  the  locations  where  the 
specified  lateral  boundaries  for  the  nest  are  in  effect. 


available.  The  nested  grids  are  rectangular  and  are  aligned 
with  the  parent  (coarser)  grid  within  which  they  are  nested. 
Additionally,  the  nested  grids  allow  any  integer  spatial  and 
temporal  refinement  ratios  of  the  parent  grid  (the  spatial  and 
temporal  refinement  ratios  are  usually  3,  but  are  not  neces¬ 
sarily  the  same  for  each  nest).  This  nesting  implementation  is, 
in  many  ways,  similar  to  the  implementations  in  the  mesos- 
cale  and  microscale  models  in  the  troposphere  and  strato¬ 
sphere;  it  allows  nesting  with  the  embedded  nested  domains 
in  the  horizontal  as  well  as  in  the  vertical  directions  (figure  2), 
and  it  allows  a  robust  novel  method  to  control  errors  near  the 
boundaries  (implicit  relaxation). 

Nested  grid  simulations  can  be  produced  using  either  1- 
way  nesting  or  2- way  nesting.  The  1-way  and  2- way  nesting 
options  refer  to  how  the  coarse  grid  and  the  fine  grid  interact. 
In  both  the  1-way  and  2-way  simulation  modes,  the  fine  grid 
boundary  conditions  (i.e.,  the  lateral  and  vertical  boundaries) 
are  interpolated  from  the  coarse  grid  forecast.  In  a  1-way  nest, 
this  is  the  only  information  exchange  between  the  grids  (from 
the  coarse  grid  to  the  fine  grid).  In  the  2- way  nest  integration, 
the  fine  grid  solution  replaces  the  coarse  grid  solution  with 
coarse  grid  points  that  lie  inside  the  fine  grid.  This  informa¬ 
tion  exchange  between  the  grids  is  now  in  both  directions 
(coarse-to-fine  for  the  fine-grid  lateral,  bottom  and  top 
boundary  computation  and  fine-to-coarse  during  the  feedback 
at  each  coarse  grid  time  step).  The  1-way  nest  setup  may  be 
run  by  one  of  two  different  methods.  One  option  is  to  produce 
the  nested  simulation  as  two  separate  simulations.  In  this 
mode,  the  coarse  grid  is  integrated  first,  and  the  coarse  grid 
forecast  is  completed.  Output  from  the  coarse  grid  integration 
is  then  processed  in  order  to  provide  boundary  conditions  for 
the  nested  run  (usually  at  a  much  lower  temporal  frequency 
than  the  coarse  grid  time  step);  and  this  is  followed  by  the 


complete  time  integration  of  the  fine  (nested)  grid.  Hence,  this 
1-way  option  is  equivalent  to  running  two  separate  simula¬ 
tions  with  a  processing  step  in  between.  The  second  1-way 
option  (lockstep  with  no  feedback),  is  run  as  a  traditional 
simulation  with  two  (or  more)  grids  integrating  concurrently, 
except  with  the  feedback  runtime  option  shut  off.  This  option 
provides  lateral  boundary  conditions  to  the  fine  grid  at  each 
coarse  grid  time  step,  which  is  an  advantage  of  the  concurrent 
1-way  method  (no  feedback). 

The  model  allows  the  refinement  of  a  coarse-grid  simu¬ 
lation  along  with  the  introduction  of  a  nested  grid.  An  option 
for  initializing  the  fine  grid  is  provided  when  using  concurrent 
1-way  and  2- way  nesting.  All  of  the  fine  grid  variables  are 
interpolated  from  the  coarse  grid.  This  option  allows  the  fine 
grid  to  start  at  a  later  time  in  the  coarse  grid’s  simulation.  A 
simulation  involves  one  outer  grid  and  may  contain  multiple 
inner  nested  grids.  Each  nested  region  is  entirely  contained 
within  a  single  coarser  grid,  referred  to  as  the  parent  grid.  The 
finer  nested  grids  are  referred  to  as  child  grids.  Using  this 
terminology,  children  are  also  parents  when  multiple  levels  of 
nesting  are  used.  The  fine  grids  may  be  telescoped  to  any  depth 
(i.e.,  a  parent  grid  may  contain  one  or  more  child  grids,  each  of 
which,  in  turn,  may  successively  contain  one  or  more  child 
grids),  and  several  fine  grids  may  share  the  same  parent  at  the 
same  level  of  nesting.  For  both  1-way  and  2- way  nested  grid 
simulations,  the  ratio  of  the  parent  horizontal  grid  distance  to 
the  child  horizontal  grid  distance  (the  spatial  refinement  ratio) 
must  be  an  integer.  In  the  cases  of  2- way  and  concurrent  1-way 
nesting,  this  is  also  true  for  the  time  steps  (the  temporal 
refinement  ratio).  The  model  does  allow  the  time  step  refine¬ 
ment  ratio  to  differ  from  the  spatial  refinement  ratio. 

The  model  uses  Arakawa-C  grid  staggering.  As  shown  in 
figure  1,  the  u,  v  and  w  components  of  the  horizontal  velocity 
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are  normal  to  the  respective  faces  of  the  grid  cell,  and  the 
thermodynamic,  density  and  chemistry  variables  are  located 
in  the  center  of  the  cell.  The  variable  staggering  has  an 
additional  column  of  u  in  the  v-direction,  an  additional  row  of 
V  in  the  y-direction  and  an  additional  row  of  w  in  the  z- 
direction  because  the  normal  velocity  points  define  the  grid 
boundaries.  The  horizontal  momentum  components  reflect  an 
average  across  each  cell  face,  while  each  thermodynamic, 
density  and  chemistry  variable  is  the  representative  mean 
value  throughout  the  cell.  Feedback  is  handled  to  preserve 
these  mean  values:  the  thermodynamic,  density  and  chemistry 
fields  are  fed  back  with  an  average  from  within  the  entire 
coarse  grid  point,  and  the  horizontal  momentum  variables  are 
averaged  along  their  respective  normal  coarse  grid  cell  faces. 
The  horizontal  interpolation  (to  instantiate  a  grid  and  to 
provide  time-dependent  lateral,  bottom  and  top  boundaries) 
does  not  conserve  mass.  The  feedback  mechanism  uses  cell 
averages  (for  thermodynamic,  density  and  chemistry  quan¬ 
tities)  and  cell-face  averages  for  the  horizontal  momentum 
fields.  The  staggering  defines  the  way  that  the  fine  grid  is 
situated  on  top  of  the  coarse  grid.  For  all  odd  ratios,  there  is  a 
coincident  point  for  each  variable:  a  location  that  has  the 
coarse  grid  and  the  fine  grid  at  the  same  physical  point.  The 
location  of  this  point  depends  on  the  variable.  In  each  of  the 
coarse  grid  cells  with  an  odd  ratio,  the  middle  fine  grid  cell  is 
the  coincident  point  with  the  coarse  grid  for  all  of  the  mass- 
staggered  fields  (figure  2).  For  the  horizontal  momentum 
variables,  the  normal  velocity  has  coincident  points  along  the 
grid  boundaries  for  odd  ratios. 

2.4.  Implicit  relaxation 

One  of  the  main  difficulties  faced  in  atmospheric  as  well  as 
oceanic  downscaling  modeling  is  the  specification  of  the 
lateral  boundary  conditions.  Usually,  the  prognostic  fields  at 
the  lateral  boundaries  of  the  nested  grid  are  specified  from  the 
large  domain.  These  fields  have  coarse  resolution,  and  are 
interpolated  in  space  and  time  to  the  nested  grid.  The 
inconsistencies  between  the  limited  and  the  large  domain 
solutions  create  spurious  reflections  that  may  propagate  and 
affect  the  solution  in  the  interior  of  the  nested  domain.  Several 
approaches  are  used  to  handle  the  lateral  boundary  conditions. 
The  flow  relaxation  scheme  ([10])  is  the  most  frequently  used 
for  atmospheric  mesoscale  forecasting  models  over  a  limited 
domain.  Lateral  open  boundary  conditions  are  often  used  in 
limited-area  ocean  modeling.  These  conditions  include  the 
radiation  condition,  the  combined  radiation  and  prescribed 
condition,  which  depends  on  the  inflow  and  outflow  regime  at 
the  boundary,  and  a  scale  selective  approach.  A  review  of 
these  methods  is  given  in  [74]. 

We  have  developed  a  novel  technique  of  implicit 
relaxation  for  limited-area  horizontally  and  vertically  nested 
ionospheric  models,  with  a  particular  emphasis  on  improved 
vertical  resolution  for  altitudes  in  the  range  of  80  km-500  km. 
The  relaxation  method  consists  of  progressively  relaxing  the 
fine  grid  fields  toward  the  coarse  grid  fields.  This  method  is 
implemented  using  an  operator  splitting  method  that  is 
applied  in  the  acoustic  steps,  in  which  a  prognostic  variable  is 
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Figure  3.  Specified  and  relaxation  zones  for  a  grid  with  a  single 
specified  row  and  column  and  four  rows  and  columns  for  the 
relaxation  zone. 

updated  first,  without  any  relaxation.  The  relaxation  is  applied 
as  an  implicit  correction  step  after  each  acoustic  step.  The 
corrected  variable  is  then  used  to  update  the  other  prognostic 
variables  ([58-60,  64,  66,  99,  100]). 

For  the  fine  grid  with  2-way  nesting  or  1-way  nesting 
(using  a  concurrent  simulation,  see  figure  2),  the  boundary 
conditions  are  specified  by  the  parent  grid  at  every  coarse  grid 
time  step.  Prognostic  variables  are  entirely  specified  in  the 
outer  row  and  in  the  column  of  the  fine  grid  through  spatial 
and  temporal  interpolation  from  the  coarse  grid  (the  coarse 
grid  is  stepped  forward  in  time  prior  to  the  advancement  of 
any  child  grid  of  that  parent).  The  nested  boundary  condition 
is  also  referred  to  as  a  relaxation,  or  nudging,  boundary 
condition.  The  nested  grid’s  lateral  and  vertical  boundaries 
are  comprised  of  both  a  specified  zone  and  a  relaxation  zone, 
as  shown  in  figure  3.  For  the  coarse  grid,  the  specified  zone  is 
determined  entirely  by  temporal  interpolation  from  the  coarse 
grid.  The  last  row  and  column  along  the  outer  edge  of  the 
nested  grid  is  entirely  specified  by  temporal  interpolation, 
using  data  from  the  coarse  grid.  The  second  region  of  the 
lateral  boundary  for  the  coarse  grid  is  the  relaxation  zone.  The 
relaxation  zone  is  where  the  model  is  nudged  or  relaxed 
towards  the  larger-scale  field  (e.g.,  rows  and  columns  2 
through  5  in  figure  3).  The  size  of  the  relaxation  zone  can  be 
adjusted  to  any  number. 

The  normal  velocity  component  that  is  located  at  the 
boundary  (u  point,  v  point  and  w  point)  is  treated  differently 
than  the  tangential  velocities  and  the  thermodynamic  vari¬ 
ables,  which  are  located  one-half  grid  point  inside  the  domain 
and  are  adjacent  to  the  boundary.  The  relaxation  boundary 
scheme  involves  smoothly  constraining  the  main  prognostic 
variables  of  the  nested  model  to  match  the  corresponding 
values  from  the  coarse  grid  model  in  a  buffer  zone  next  to  the 
boundary  called  the  ‘relaxation’  zone.  The  flow  relaxation 
scheme  is  a  combination  of  Newtonian  and  diffusive 
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relaxation,  and  has  the  form: 

d,w  =-N{x){\if-  y/")  +  D(x)d^{y/-  y/") 

where  i/a  is  a  prognostic  variable  of  the  limited-area  model 
that  needs  to  be  relaxed  to  the  corresponding  variable  from 
the  coarse  grid  model  y/^  x  denotes  the  direction  normal  to 
the  boundary  and  N(x)  and  D(x)  are  the  Newtonian  and  dif¬ 
fusive  relaxation  factors.  The  choice  of  the  profiles  and  the 
values  of  the  coefficients  N(x)  and  D(x)  in  the  relaxation 
zones  control  the  reflections  at  the  boundary. 

The  relaxation  is  implemented  after  each  acoustic  step  as 
an  implicit  correction.  We  allow  to  denote  the  pertur¬ 

bation  of  the  updated  variable  after  each  acoustic  time  step. 
The  relaxation  is  then  applied  as  a  correction  in  a  subsequent 
step  using  the  following  implicit  flow  relaxation  equation: 


-A^(i 

1 

+ 

{(cr+vL-w;::') 

-l{y 

-  <■") 

+  (v/ 

where  5t  is  the  acoustic  time  step,  5x  is  the  grid  spacing  and 
is  the  coarse  grid  value  that  is  interpolated  in  space  and 
to  the  time  step  (n+1).  i/a"  '"'^^  =  i/a/''^^  —  i/a/,  where  is  the 
total  fine  grid  value,  i/a^  is  the  most  updated  value  in  the  RK 
step  and  i/a"  '''^^  is  the  perturbation,  with  respect  to  i/a/.  For  the 
prognostic  variables,  located  at  half-grid  points  that  are 
adjacent  to  the  lateral  boundary,  this  implicit  equation  is 
solved  for  i/a"  '"'*'^  along  the  relaxation  zone,  subject  to  the 
boundary  conditions: 


(a)  (b) 


Figure  4.  Vertical  cross  sections  of  density  irregularities  (color)  from 
large  domain  simulations  at  three  different  times:  (a)  tl  =  191,  (b) 
t2  =  291  and  (c)  t3  =  491.  The  two  vertical  dashed  lines  enclose  the 
domain  that  is  used  for  limited-area  simulations. 


.,n+l  ^ 


...1  ^  cn.l  _  . 


where  5’  is  the  index  of  the  last  relaxed  point  in  the  interior  of 
the  domain.  For  the  normal  velocities  at  the  boundary,  the 
same  system  will  be  solved,  except  that  consistency  between 
the  coarser  and  the  finer  mass  fluxes  in  the  continuity 
equation  will  be  imposed;  that  is: 


V  •  y  =  y  •  y' 


where  y  and  V""  are  the  velocity  vectors  from  the  nested  and 
the  coarse  grids.  Since  the  tangential  velocities  that  are 
adjacent  to  the  boundary  are  imposed  by  the  coarse  grid 
model,  the  above  relation  reduces  to: 

dU  _  dW 
dx  dx 


where  U  and  V  are  the  normal  components  of  y  and  v\ 
respectively.  This  relationship  is  imposed  implicitly  at  the 
lateral  boundaries,  and  the  implicit  equation  is  solved  for 
along  the  relaxation  zone  as  stated  above,  but  with  the 


following  implicit  boundary  conditions: 

jj/f,n+l  _  ^jf/,n+l 

^5+1  ~  ^  5+1 


and 


U«.n+l  _  ^ 

The  Newtonian  and  diffusive  relaxation  times  are  fixed 
by  the  choice  of  the  coefficients  and  T)  .  Optimal  profiles 
for  the  Newtonian  relaxation  coefficient  N  are  proposed  in 
[10]  and  [55].  They  are  constructed  in  such  a  way  that,  under 
idealized  conditions,  the  unwanted  partial  reflection  of  the 
outgoing  waves  that  leave  the  domain  is  minimized.  More 
recently  [71],  conducted  a  detailed  theoretical  and  numerical 
study  on  the  choice  of  the  Newtonian  and  diffusive  relaxation 
coefficients,  in  which  they  compared  different  profiles  and 
relaxation  times.  They  concluded  that  a  diffusive  relaxation, 
combined  with  exponential  or  optimized  profiles,  gives  better 
results.  Also,  they  presented  a  formula  for  what  they  called  the 
leading  coefficients  and  which  they  derived  using  the 
criterion  of  minimum  reflection.  These  are  nondimensional 
coefficients  at  the  first  relaxed  grid  point.  These  coefficients. 
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together  with  the  relaxation  profiles,  determine  the  entire 
relaxation  time  in  the  relaxation  zone  through  the  relations: 


N.  =  — —N*N.  and  D.  =  — —D*D: 
‘  2Ax  ^  ‘  ‘  2Ax  ^  ' 


In  these  relations,  c  is  the  phase  speed  of  the  fastest  wave 
and  N-  and  D.  are  the  normalized  profiles  of  the  coefficient 
(N^  =  1  and  =  1).  In  our  implicit  relaxation,  we  choose  a 
five-point  or  nine-point  deep  relaxation  zone,  in  which  the 
Newtonian  and  diffusive  relaxation  times  are  applied.  Fol¬ 
lowing  [71],  we  choose  =  0.9  and  =  0.9;  c  is  the 
speed  of  sound,  and  an  optimized  profile  is  computed  for  N. 
and  D.  using  the  algorithm  presented  in  [55].  These  specifi¬ 
cations  entirely  determine  the  Newtonian  and  diffusive 
relaxation  times. 

Figure  4  shows  a  test  of  the  implicit  relaxation  for  an 
idealized  simulation.  It  shows  cross-sections  of  density  irre¬ 
gularities  from  a  large  domain  simulation,  using  open 
boundary  conditions.  The  grid  spacing  is  1.5  km.  The  cross- 
sections  are  shown  for  three  different  times  tl  =  191,t2  =  291 
and  t3  =  491,  after  density  irregularities  have  developed.  The 
two  vertical  dashed  lines  enclose  the  nested  domain,  which  is 
used  for  limited-area  simulations.  The  nested  domain  is 
chosen  to  be  outside  the  main  density  anomaly  region.  This 
simulation  uses  a  grid  spacing  that  is  three  times  finer  than  the 
one  used  for  the  large  domain  simulation.  All  the  fields, 
including  the  densities,  are  relaxed  to  equal  the  fields  from  the 
large  domain  by  using  5-point-wide  implicit  relaxation.  These 
tests  show  that  the  density  irregularity  and  the  associated 
dynamics  smoothly  cross  the  boundary  without  noticeable 
reflections.  The  Newtonian  and  diffusive  relaxation  schemes 
used  in  the  implicit  scheme  are  fixed  by  the  choice  of  the 
optimal  coefficients  that  are  constructed  in  such  a  way  that  the 
unwanted  partial  reflections  of  the  outgoing  waves  that  leave 
the  domain  is  minimized. 

We  also  mention  resolving  techniques  for  front  tracking 
and  density  current  fronts  that  have  been  used  in  the  context 
of  neutral  fluids  in  order  to  study  classical  Rayleigh-Taylor 
instabilities  ([25])  and  the  dynamics  of  fronts  of  density 
currents  that  are  associated  with  geophysical  flows  ([35,  67]). 
These  efficient  computational  methods,  previously  used  for 
neutral  fluids,  are  extended  to  charged  interfaces/fronts  in  the 
context  of  ionospheric  dynamics. 


3.  Multiscale  nested  simulations  of  primary  and 
secondary  Rayleigh-Taylor  instabilities  in 
ionospheric  flows,  equatorial  spread  F  (ESF) 

Plasma  irregularities  and  inhomogeneities  in  the  F  region, 
caused  by  plasma  instabilities,  manifest  as  spread  F  echoes. 
The  scale  sizes  of  the  density  irregularities  range  from  a  few 
meters  to  a  few  hundred  kilometers,  and  the  irregularities  can 
appear  at  all  latitudes.  However,  spread  F  echoes  in  the 
equatorial  region  can  be  particularly  severe.  At  night,  fully 
developed  spread  F  is  characterized  by  plasma  bubbles,  which 
are  vertically  elongated  wedges  of  depleted  plasma  that  drift 
upward  from  beneath  the  bottomside  F  layer.  Under  certain 


conditions,  a  density  perturbation  can  trigger  the  Rayleigh- 
Taylor  instability  on  the  bottomside  of  the  F  layer.  Once 
triggered,  density  irregularitues  develop,  and  the  field-aligned 
depletions  then  bubble  up  through  the  F  layer.  The  east-west 
extent  of  a  disturbed  region  can  extend  up  to  several  thousand 
kilometers,  with  the  horizontal  distance  between  the  separate 
depleted  regions  extending  to  tens  to  hundreds  of  kilometers. 
The  plasma  density  in  the  bubbles  can  be  up  to  two  orders  of 
magnitude  lower  than  that  in  the  surrounding  medium.  When 
spread  F  ends,  the  upward  drift  ceases  and  the  bubbles 
become  fossilized. 

The  evolution  of  the  ESF  is  a  strongly  nonlinear  phe¬ 
nomenon,  with  multiscale  interactions  for  ionospheric 
dynamics.  The  large-scale  primary  Rayleigh-Taylor  mode 
can  promote  a  hierarchy  of  smaller-scale  plasma  instability 
processes  that  give  rise  to  a  wide  spectrum  of  irregularities. 
The  presence  of  these  small  scale  irregularities  was  evidenced 
from  the  observations  that  showed  the  coexistence  of  kilo¬ 
meter  disturbances  and  meter  scale  disturbances  in  the 
nighttime  equatorial  F  region  ([7]).  Thus,  simulations  that  use 
coarse  and  fixed  resolutions  cannot  resolve  the  small  scale 
disturbances.  Poor  resolution  of  these  scales  can,  in  turn, 
affect  the  accuracy  of  the  larger  scales  due  to  nonlinearity. 
Obviously,  one  can  design  a  computer  model  with  a  very  high 
spatial  resolution  everywhere.  But  then,  the  simulations  will 
be  prohibitively  expensive,  particularly  in  three  dimensions. 

Considerable  efforts  have  been  made  in  the  development 
of  computer  models  that  simulate  the  generation  and  evolution 
of  Equatorial  Spread  F  (ESF)  dynamics.  Most  of  these  models 
use  periodic  boundary  conditions  with  fixed  spatial  resolutions 
in  both  the  horizontal  and  the  vertical  directions.  The  transport 
algorithm  that  is  commonly  used  in  these  models  is  the  flux- 
corrected  transport  method  ([120]).  This  scheme,  however, 
tends  to  develop  relatively  high  numerical  diffusivity  that  can 
be  comparable  to  the  physical  diffusivity  produced  by  the 
subgrid  scale  turbulent  mixing  [27].  The  Crank-Nicholson  (C- 
N)  scheme  was  also  used  in  some  ionospheric  models  ([4,  32]). 
The  C-N  scheme,  however,  can  produce  excessive  noise, 
thereby  degrading  the  numerical  solution. 

In  this  section,  nested  numerical  simulations  of  iono¬ 
spheric  plasma  density  structures  that  are  associated  with  the 
nonlinear  evolution  of  the  Rayleigh-Taylor  (RT)  instability  in 
the  Equatorial  Spread  F  (ESF)  are  presented.  The  numerical 
implementation  of  the  nested  model  uses  a  spatial  dis¬ 
cretization  with  a  C  grid  staggering  configuration,  in  which 
normal  velocities  of  ions  and  electrons  are  staggered  one-half 
grid  length  from  the  density  of  the  charged  particles  (figures  1 
and  2).  The  advection  of  charged  particles  is  computed  with  a 
fifth-order,  accurate  in  space.  Weighted  Essentially  Non- 
Oscillatory  (WENO)  scheme.  The  continuity  equation  is 
integrated  using  a  third-order  Runge-Kutta  (RK)  time  inte¬ 
gration  scheme.  The  equation  for  the  electric  potential  is  solved 
at  each  time  step  with  a  multigrid  method.  For  the  limited-area 
and  nested  simulations,  the  lateral  boundary  conditions  are 
treated  via  implicit  relaxation,  which  is  applied  in  the  buffer 
zones  where  the  density  of  the  charged  particles  for  each  nest  is 
relaxed  to  the  density  of  the  charged  particles  obtained  from  the 
parent  domain  (figure  3).  High  resolution  in  targeted  regions. 
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offered  by  the  nested  model,  is  able  to  resolve  secondary  RT 
instabilities  and  to  improve  the  resolution  of  the  primary  RT 
bubble  compared  to  the  coarser  large  domain  model.  The 
computational  results  are  validated  via  a  large  domain  simu¬ 
lation,  in  which  the  resolution  is  increased  everywhere. 

The  equatorial  ionosphere  is  considered  in  a  slab  geo¬ 
metry,  with  Earth’s  magnetic  field  B  directed  along  the  y  axis. 
We  choose  the  unit  vectors  i,  j  =  B/B  and  k  to  be  directed 


along  the  westward,  northward  and  upward  directions, 
respectively.  A  set  of  plasma  two  fluid  equations,  which 
describe  the  conservation  of  momentum,  mass  and  current, 
are  considered.  These  equations  govern  the  motion  of  the 
species  5'  ions  (s  =  i)  and  the  electrons  (s  =  e),  the  masses  of 
which  are  represented  as  Mi  and  Me,  respectively.  The 
computational  methodology  that  is  presented  in  section  2  is 
validated  via  a  large  domain  simulation  of  the  interaction  of 


Longitude  (km) 


_i _ i _ 1- 

-so  0  so 

Longitude  (km) 


I 


a; 

-U 

:3 


< 


(a) 


(b) 


so  200 


latitude  (km) 


(c) 


(d) 


Figure  5.  3D  large  domain  simulation  of  the  ESF.  The  neutral  wind  is  125  m  s“\  and  the  imposed  external  electric  field  is  given  by  EJB=  100, 
where  B  is  the  magnetic  field.  The  horizontal  axis  represents  the  east-west  range  (a)  and  (b),  and  the  north-south  range  (c)  and  (d).  The  solid 
curves  represent  the  iso-density  contours  after  2000  s.  The  simulation  is  initialized  with  a  small  3D  density  perturbation  that  is  superimposed  with 
the  background  density  profile.  The  x-z  cross-sections  are  at  (a)  y  =  0km  and  (b)  y=  109  km.  The  y-z  cross-sections  are  at  (c)  at  x  =  0  and 
(d)  X  =  3 1  km. 
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Figure  6.  (a)  Large  domain  simulation  of  the  ESF.  The  neutral  wind  magnitude  is  125  m  s“\  and  the  imposed  external  electric  field  is  given 
by  EJB  =  100,  where  B  is  the  magnetic  field.  The  horizontal  axis  represents  the  east-west  range.  The  cross-section  is  taken  at  y  =  0.  The  solid 
curves  represents  the  iso-density  contours  after  2000  s.  The  simulation  is  initialized  with  a  small  density  perturbation  that  is  superimposed  on 
the  background  density  profile.  The  two  dashed  vertical  lines  represent  the  horizontal  boundary  of  the  domain  used  in  the  limited  area  and  the 
nested  simulations,  (b)  Limited  area  domain  simulation  of  the  ESF.  The  lateral  boundary  conditions  use  the  implicit  relaxation  that  is  applied 
in  buffer  zones  where  the  density  of  the  charged  particles  is  relaxed  to  equal  the  density  of  the  charged  particles  obtained  from  the  parent 
domain,  shown  in  (a).  The  horizontal  axis  represents  the  east- west  range.  The  solid  curves  represent  the  iso-density  contours  after  2000  s. 
The  simulation  is  initialized  with  a  small  density  perturbation  that  is  superimposed  with  the  background  density  profile.  The  horizontal  grid 


spacing  is  the  same  as  the  spacing  in  the  simulation  shown  in  (a). 

the  equatorial  ionosphere  with  the  eastward  F  region  neutral 
wind  in  the  presence  of  evolving  spread  F  bubbles.  Standard 
background  profiles  of  ion-neutral  collision  frequency, 
recombination  and  the  background  electron  density  profile  as 
a  function  of  altitude  are  used  in  the  numerical  simulations 
([58,  99,  100]).  An  assessment  of  the  performance  of  local 
ensemble  Kalman  filter  data  assimilation  techniques  applied 
to  the  ionospheric  flows  is  given  in  [16]. 

The  extent  of  the  large  domain  simulation  is  [-200  km, 
200  km]  in  the  horizontal  direction  and  [300  km,  550  km]  in 
the  vertical  direction,  with  a  grid  spacing  of  4  km,  and  2.5  km 
in  horizontal  directions  and  vertical  directions,  respectively. 
Figure  5  shows  a  3D  large  domain  simulation  of  the  Equa¬ 
torial  Spread  F  (ESF).  The  neutral  wind  is  125  m  s“\  and  the 
imposed  external  electric  field  is  given  by  Eo/B  =  100,  where 
B  is  the  magnetic  field.  The  horizontal  axis  represents  the 
east-west  range,  (a)  and  (b),  and  the  north-south  range,  (c)  and 
(d).  The  solid  curves  represent  the  iso-density  contours  after 
2000  s.  The  simulation  is  initialized  with  a  small  3D  density 
perturbation  superimposed  to  the  background  density  profile. 
The  x-z  cross-sections  are  at  (a)  y  =  0  km  and  (b)  y  =  109  km. 
The  y-z  cross-sections  are  at  (c)  at  x  =  0  and  (d)  x  =  3 1  km. 
The  development  of  the  ESF  as  a  large-scale  bubble  is  evi¬ 
dent.  The  top  of  this  bubble  reaches  a  high  altitude  and  is 
located  below  500  km.  We  note  that  our  model  does  not 


develop  spurious  small-scale  numerical  noise,  which  occurred 
in  the  Crank-Nicholson  scheme  in  [5]. 

To  test  the  performance  of  the  lateral  implicit  relaxation, 
the  outputs  from  the  large  domain  are  truncated  within  a 
limited  area  [-125  km,  150  km],  which  is  represented  by  the 
dashed  vertical  lines  in  figure  6;  and  a  second  simulation  is 
conducted  using  the  same  model  within  this  limited  area 
without  increasing  the  resolution  (with  4  km  horizontal  grid 
spacing).  This  simulation  uses  the  same  configuration  that 
was  used  in  the  large  domain  (same  horizontal  and  vertical 
grid  distributions),  except  that  the  number  of  grid  points  in  the 
horizontal  is  reduced;  and  both  the  electric  potential  and 
plasma  density  are  specified  at  the  lateral  boundaries.  The 
density  computed  from  the  nested  model  is  relaxed  to  the 
coarse  outputs  from  the  large  domain  by  using  the  implicit 
relaxation.  The  width  of  the  relaxation  zone  is  five  points. 
Figure  6  shows  the  field  of  the  iso-density  contours,  simulated 
by  the  limited-area  model  after  2000  s.  The  simulation  results 
from  the  large  domain  (figure  5)  and  from  the  limited-area 
domain  (figure  6)  are  indistinguishable  over  the  limited-area 
domain,  thereby  demonstrating  that  the  nested  model  and  the 
method  proposed  for  the  relaxation  perform  well. 

To  test  the  impact  of  the  resolution  on  the  evolution  of  the 
ESF,  a  nested  simulation  was  conducted  using  the  same  setup 
as  the  one  used  for  the  limited  area  simulation,  but  the 
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Figure  7.  (a)  Nested  domain  simulation  of  ESF.  The  resolution  is  doubled  (in  the  finest  nest  only).  The  lateral  boundary  conditions  use  the 
implicit  relaxation  that  is  applied  in  the  buffer  zones  where  the  density  of  the  charged  particles  is  relaxed  to  equal  the  density  of  the  charged 
particles  interpolated  from  the  parent  domain.  The  solid  curves  represent  the  iso-density  contours  after  2000  s.  The  simulation  is  initialized 
with  a  small  density  perturbation  that  is  superimposed  on  the  background  density  profile.  Note  that  the  nested  simulation  resolves  the 
secondary  RT  instabilities.  Also,  the  primary  disturbance  is  more  developed  and  reaches  higher  altitudes,  (b)  Large  domain  simulation  with 
high  resolution  (the  resolution  is  doubled  everywhere).  This  demonstrates  that  the  secondary  RT  instabilities  (bubbles)  that  are  resolved  by 
the  nested  simulations  are  not  artifacts,  and  that  the  nesting  offers  more  accurate  results  at  low  computational  cost. 


resolution  is  doubled  (the  horizontal  grid  spacing  is  2  km).  The 
lateral  boundary  conditions  use  the  imphcit  relaxation  that  is 
applied  in  the  buffer  zones  where  the  density  of  the  charged 
particles  of  the  nest  is  relaxed  to  equal  the  density  of  the 
charged  particles  of  the  nest  that  are  interpolated  in  time  and 
space  from  the  parent  domain,  as  described  in  section  2. 
Figure  7  shows  the  field  of  iso-density  contours,  simulated  by 
the  nested  model  after  2000  s.  The  results  show  that  in  addition 
to  the  main  spread  F  bubble  that  was  also  resolved  in  the  large 
domain  simulation,  there  is  a  generation  of  secondary  Ray- 
leigh-Taylor  instabilities  or  secondary  bubbles  that  were  not 
resolved  in  the  coarse  grid  simulation.  In  addition,  the  primaiy 
(the  large  scale)  disturbance  is  more  developed  and  reaches 
higher  altitudes. 

To  show  that  the  secondary  ESF  instability  and  the 
higher  altitude  of  the  top  of  the  main  bubble  are  not  artifacts 
due  to  numerical  errors  in  the  nested  model,  another  simu¬ 
lation  was  conducted. 

We  used  a  large  domain,  as  was  also  used  in  the  first 
simulation;  however,  in  the  second  simulation,  the  resolution  was 
doubled  everywhere.  The  results  of  this  simulation  are  presented 
in  figure  8.  We  noted  that  by  increasing  the  resolution  of  the 
parent  model,  the  secondary  instability  was  now  resolved,  as  in 
the  nested  simulation  (figure  7).  Also,  the  primaiy  disturbance  is 
more  developed  and  reaches  a  higher  altitude  (above  500  km)  as 
the  resolution  is  doubled.  Thus,  the  high  resolution  in  the  targeted 


regions  offered  by  the  nested  model  is  found  to  be  critical  for  the 
resolution  of  the  small-scale  ionospheric  plasma  density  struc¬ 
tures  and  for  the  improvement  of  the  large-scale  bubble,  which 
are  associated  with  the  evolution  of  the  primary  and  secondary 
Rayleigh-Taylor  instabihties  in  the  Equatorial  Spread  F. 


4.  Impacts  of  time  varying  neutrai  winds  on  fine 
scaie  structure  of  sporadic  E  and  intermediate 
iayers 

In  this  section,  we  present  three-dimensional  numerical  stu¬ 
dies  of  the  E  and  lower  F  region  ionosphere,  coupled  with  the 
neutral  atmosphere  dynamics.  The  inclusion  of  neutral 
dynamics  in  the  model  allows  us  to  study  the  transport  pat¬ 
terns  of  plasma  density  and  to  examine  the  charge  of  the 
neutral  interactions  over  the  full  evolution  cycle  of  an  inertial 
gravity  wave  when  the  background  flow  spins  up  from  rest, 
saturates  and  eventually  breaks.  Using  Lagrangian  analyses, 
we  show  the  mixing  patterns  of  the  ionospheric  responses, 
formation  and  subsequent  evolution  of  nonequilibrium  iono¬ 
spheric  layers.  The  corresponding  plasma  density  in  this  flow 
develops  complex  wave  structures  and  small-scale  patches 
during  the  gravity  wave  breaking  event. 

The  lower  ionospheric  altitudes  are  challenging  to  model 
because  they  are  too  low  for  orbiters  and  too  high  for 
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Figure  8.  Nesting  summary:  high  resolution  modeling  of  the  multiscale  ionospheric  dynamics  over  a  limited  spatial  range  in  the  altitude, 
latitude  and  longitude,  which  provides  accurate  results  at  a  low  computational  cost,  (a)  With  nesting  and  implicit  relaxation  to  the  coarse 
model  (results  verified  by  high  resolution  computations  on  a  large  domain)  and  (b)  without  nesting. 


radiosondes  to  take  direct  measurements.  In  recent  years, 
computer  simulations  of  the  Earth’s  ionosphere  have  become 
a  prevailing  tool  used  to  obtain  properties  of  plasma  flows  in 
the  ionosphere,  especially  at  low  altitudes.  Numerical  models 
have  been  developed  to  study  Rayleigh-Taylor  instabilities  in 
the  equatorial  spread  F  and  ionospheric  responses  to  neutral 
atmospheric  motions  in  mid-latitude.  In  these  studies,  the 
neutral  dynamics  are  prescribed  as  idealized  velocity  fields, 
such  as  a  constant  drift  flow  or  empirical  shear  flow  models; 
and  they  are  specifically  used  for  mid-latitude  simulations, 
linear  models  of  inertial  gravity  waves  (IGW)  or  data  sets 
from  other  atmospheric  models,  which  are  influenced  by 
many  dynamical  processes. 

We  study  3D  dynamics  in  the  E-  and  F-region  of  the 
ionosphere  due  to  interactions  with  neutral  winds.  A  paral¬ 
lelized  DNS  solver  is  used  to  solve  the  coupled  system  of 
charged  flows  and  the  3D  Navier-Stokes  equations  for  neu¬ 
trals.  In  the  current  literature,  the  background  neutral  flow  is 
usually  prescribed  and  not  coupled  into  the  dynamics.  We 
include  this  coupling  permitting  time-dependent  neutral 
winds  in  the  model  ([58,  99,  100]).  Multi-nested  limited-area 
high  resolution  simulations  of  the  ionosphere  for  a  range  of 
altitudes  from  80  km  to  500  km,  which  resolve  the  coupling 
between  time-dependent  neutral  winds  and  ions/electrons  at 
mesoscale/microscale  resolutions,  are  conducted. 

Sharply  defined  layers/density  interfaces  determine  fine 
scale  structure  of  the  ionosphere.  Charged  shear  layers  form 
along  contours  that  are  defined  by  very  large  gradients  in 
densities  and  scintillations.  Inhomogeneous  turbulent 


dynamics,  which  are  adjacent  to  these  layers,  stretch  and 
distort  the  larger  scales  in  order  to  counter  the  tendency  of  the 
interface  to  thicken.  Our  focus  is  on  the  multiscale  dynamics 
of  the  flow  within  and  near  the  density  interfaces  (electron 
density/ionization  fronts).  The  root  mean  square  values  of 
density  and  velocity,  and  also  their  flatness  factors,  are  given 
by  statistical  averaging  across  these  layers,  taking  into 
account  the  degree  to  which  the  thin  layers  are  sufficiently 
convoluted  to  be  ‘space-filling’,  with  the  thickness  of  the 
layer  edges  being  approximately  equal  to  the  Taylor  micro¬ 
scale.  Through  the  blocking  of  the  external  scale  eddies, 
which  occur  outside  the  layers  as  they  impinge  onto  the 
interfaces,  upscale  and  downscale  energy  transfer  processes 
lead  to  distorted  inertial-range  motions  with  a  wide  range  of 
length  scales  and  characteristic  self-similar  power-law  spec¬ 
tra.  Here,  the  skewness  of  the  velocity  derivatives  becomes 
negative.  The  local  dynamics  near  such  density  interfaces 
(inhomogeneous,  non-Kolmogorov)  is  quite  different  to 
dynamics  in  which  a  wide  spectrum  of  eddies  interact 
homogeneously  through  the  flow  (homogeneous,  isotropic 
Kolmogorov  turbulence).  The  structure  of  these  thin  layers  is 
identified  by  conditionally  averaging  the  fields  relative  to  the 
moving  interface.  The  shear  layer  itself  is  unstable  and  gen¬ 
erates  eddies  on  the  scale  of  the  shear  layer  thickness.  These 
strongly  influence  the  average  movement  of  the  interface. 
Simulations  show  how  the  interface  moves  outward  at  the 
boundary  entrainment  velocity.  The  very  high  velocity  gra¬ 
dients  across  the  thin  layer  are  maintained  by  the  interaction 
between  the  large-scale  impinging  eddies  and  the  large-scale 
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shear  outside  the  interface.  Note  that  conventional  turbulence 
models,  even  when  used  for  unsteady  flows,  simulate  the 
energetic  diffusive  action  of  small-scale  eddies  and  not  the 
counter-diffusion  action  of  highly  inhomogeneous  larger- 
scale  eddies  near  the  interfaces.  Even  if  the  initial  third 
moments  (or  skewness)  of  the  large  and  small  scale  gradients 
are  zero,  these  interactions  generate  a  negative  third  moment 
of  the  derivatives.  The  mechanism  also  leads  to  the  third 
moment  of  the  velocity  fluctuations  being  negative  as  a  result. 
The  very  high  velocity  gradients  that  separate  regions  of  flows 
with  differing  levels  of  densities  have  a  complex  geometry. 

Polarization  charges  tend  to  accumulate  at  sharply 
defined  density  interfaces  ([99,  100]).  As  many  researchers 
have  commented,  new  concepts  are  needed  in  order  to 
understand  the  spatially  intermittent  structure  of  turbulent 
interfaces.  As  for  other  well-defined  large-scale  coherent 
structures,  a  combination  of  quasi-deterministic  modeling  of 
the  defined  types  of  the  interface  structure  and  the  probabil¬ 
istic  analysis  of  their  occurrence  can  describe  the  character¬ 
istic  features  as  well  as  their  overall  statistics.  This  approach 
can  overcome  substantial  disadvantages  of  the  current  statis¬ 
tical  modeling  approach  by  providing  physical  explanations 
of  the  key  features  of  charged  turbulent  interfaces,  namely, 
the  extreme  non-Gaussianity  of  its  smallest  scales  and  its 
critical  upscale  transfer  processes,  which  lead  to  the  formation 
of  internal  shear  layers.  High  resolution  simulations,  with 
combined  Eulerian  and  Lagrangian  techniques,  allow  us  to 
resolve  the  sharpness  of  the  charged  density  interfaces,  which 
are  often  overlooked. 

The  fine  scale  structure  of  the  3D  time-mean  electron/ion 
density  distributions  are  determined  by  time-mean  vertical  and 
horizontal  electron/ion  density  fluxes:  (//w'),  {p^u'),  (pjv'). 
Here,  electron/ion  densities  and  velocities  are  decomposed  into 
ensemble  means  and  fluctuations  p  =  {p)  +  {pj)^  u  =  {u)  +  u', 
V  =  (v)  -Hv',  w  =  {w)-\-w'.  The  effects  of  ionospheric  processes 
on  the  variance  of  spatial  density  anomalies  at  ionospheric 
meso/micro  scales  need  to  be  assessed.  If  a  process  generates 
spatial  density  anomalies,  whether  they  are  positive  or  nega¬ 
tive,  it  leads  to  an  increase  in  {p'Y.  There  are  several  processes 
capable  of  affecting  the  electron/ion  density  fields.  One  is 
related  to  electron/ion  vertical  motions.  Here,  the  question  of 
primary  interest  is  the  sign  of  the  term  {p'w'),  normalized  by 
variances  of  p'  and  w'.  Another  process  capable  of  affecting  the 
3D  density  distribution  is  related  to  density  fluxes,  which  are 
induced  by  eddies  that  are  associated  with  neutral  velocities. 
The  term  {p^w')  represents  the  generation  of  spatial  electron/ 
ion  density  anomalies  through  vertical  motions;  (p'Y  are 
generated  by  upward  motions  in  regions  with  higher  electron/ 
ion  densities  or  by  downward  motions  in  regions  with  lower 
densities.  Similarly,  the  terms  {p^u')  and  (pjv')  are  responsible 
for  horizontal  fluxes  and  are  of  particular  interest  in  assessing 
the  relative  importance  of  various  processes  used  in  the  gen¬ 
eration  and  the  destruction  of  spatial  electron/ion  density 
anomahes.  We  investigate  the  generation  and  the  destruction  of 
spatial  density  anomalies  through  the  ensemble  and  time-mean 
fluxes  associated  with  the  mean  vertical  and  horizontal  motions 
and  through  the  density  fluxes  induced  by  time  dependent 
velocities. 


We  note  that  vertical  motions  can  alter  the  electron 
density  distribution  at  ionospheric  meso/microscales.  We 
compare  the  distribution  of  {w')  with  that  of  (pj)  to  identify 
altitudes  in  which  the  lower  electron  density  regions  sink  and 
altitudes  in  which  the  higher  electron  density  regions  rise. 
These  vertical  motions  lead  to  a  generation  of  spatial  electron 
density  anomalies.  When  integrated  over  the  area  of  the 
density  layer,  a  net  generation  of  spatial  density  anomalies 
through  time-mean  vertical  motions  can  be  found.  The  ver¬ 
tical  profile  of  (pjw'),  integrated  over  each  density  layer,  is  of 
particular  interest.  We  identify  altitudes  where  large  con¬ 
tributions  come  from  the  vertical  density  flux.  Similar  to 
vertical  motions,  horizontal  density  fluxes  {p^u')  and  (pjv'), 
which  are  induced  by  u'  and  v',  can  affect  the  electron/ion 
density  fields  in  the  mesoscale/microscale  domain.  Detailed 
temporal  and  spatial  characteristics  of  fluctuating  density 
interfaces  (electron  density  and  ionization  fronts)  were 
investigated  in  [58,  99,  100]. 

Sporadic  E  layers  are  ionization  enhancements  in  the  E 
region  at  altitudes  between  90  and  120  km.  The  layers  tend  to 
occur  intermittently  and  can  be  seen  at  all  latitudes.  Sporadic 
E  layers  at  mid-latitudes  are  primarily  a  result  of  wind  shears 
([88]),  but  they  can  also  be  created  by  diurnal  and  semi¬ 
diurnal  tides.  The  layers  are  formed  when  the  vertical  ion  drift 
changes  direction  with  the  altitude,  and  the  layers  occur  at  the 
altitudes  where  the  ion  drift  converges.  In  the  E  region,  the 
zonal  neutral  wind  is  primarily  responsible  for  inducing  ver¬ 
tical  ion  drifts,  which  result  from  a  dynamo  action.  A  reversal 
of  the  zonal  neutral  wind  with  the  altitude  results  in  ion 
convergence  and  divergence  regions.  The  helical  nonparallel 
neutral  flows  are  frequently  observed  in  the  lower  ionosphere 
at  mid-latitudes.  Three-dimensional  shear  instabilities,  which 
are  attributed  to  the  breakdown  of  polarized  wind  fields,  are 
common  at  ionospheric  altitudes  that  encompass  the  E  layers. 
These  flows  are  characterized  by  neutral  wind  profiles  with 
intense  variations  in  both  eastward  and  westward  components 
over  a  short  vertical  distance.  This  intense  nonparallel  shear 
can  lead  to  the  development  of  3D  helical  Kelvin-Helmholtz 
instability,  with  characteristics  that  are  different  from  the  one 
found  in  the  parallel  flows.  The  stability  criterion  for  a  general 
nonparallel  polarized  wind  field  has  been  recently  established 
in  [56].  Contrary  to  the  parallel  flows,  this  criterion  is  ani¬ 
sotropic  and  depends  on  the  direction  of  the  wave  vectors  of 
the  unstable  modes.  The  physical  and  dynamical  processes 
associated  with  polarized  instabilities  are  inherently  three- 
dimensional.  The  horizontal  velocity  vector  (f7(v3)  and 
y(v3))  (horizontally  and  time  averaged)  rotates  with  the 
vertical  coordinate  (j^).  We  allow  a(x^)  to  denote  the  angle 
between  the  vector  horizontal 

wave  vector.  Then,  we  can  obtain  the  following  expression 
for  the  polarized  plasma  Richardson  number. 


Here,  the  definition  of  N  is  based  on  the  number  density  n 
instead  of  the  temperature,  The  polarized  plasma 
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Figure  9.  Heat  flux  shows  strong  vertical  mixing. 


Richardson  number  takes  into  account  horizontal  anisotropy,  the 
angle  between  the  horizontal  wave  vectors  and  the  velocity  vectors 
at  each  vertical  level.  In  the  case  of  polarized  wind  fields  in  stably 
stratified  neutral  environments,  the  polarized  Richardson  number 
was  rigorously  introduced  in  [56,  66].  The  above  expression 
generahzes  this  definition  to  ionospheric  flows,  in  which  the 
electron  number  density  is  used  in  place  of  the  temperature,  since 
the  ionosphere  is  stratified  with  respect  to  the  electron  density. 
KeUey  ([50],  Page  159)  discusses  a  plasma  Richardson  number,  in 
the  context  of  parallel  flows  in  the  ionosphere,  and  the  role  of 
velocity  shear  in  triggering  convective  ionospheric  storms.  The 
above  expression  generalizes  his  definition. 

In  contrast  to  sporadic  E  layers,  intermediate  layers  are 
broad  (10-20  km  wide)  and  occur  in  the  altitude  range  of 
120-180  km.  They  frequently  appear  at  night  in  the  valley 
between  the  E  and  E  regions,  but  they  can  also  appear  during 
the  day.  They  tend  to  form  on  the  bottomside  of  the  E  region 
and  then  slowly  descend  throughout  the  night  toward  the  E 
region.  As  with  sporadic  E  layers,  intermediate  layers  can 
occur  at  all  latitudes,  have  a  large  horizontal  extent,  and 

In  contrast  to  sporadic  E  layers,  intermediate  layers  are 
broad  (10-20  km  wide)  and  occur  in  the  altitude  range  of 
120-180  km.  They  frequently  appear  at  night  in  the  valley 
between  the  E  and  E  regions,  but  they  can  also  appear  during 
the  day.  They  tend  to  form  on  the  bottomside  of  the  E  region 
and  then  slowly  descend  throughout  the  night  toward  the  E 
region.  As  with  sporadic  E  layers,  intermediate  layers  can 
occur  at  all  latitudes,  have  a  large  horizontal  extent,  and  have 
an  order  of  magnitude  density  enhancement  relative  to  the 
background  densities  ([88]).  Intermediate  layers  are  primarily 
a  result  of  wind  shears  connected  with  the  semi-diurnal  tide. 
In  the  E-E  region  valley  (130-180  km),  the  meridional  neutral 
wind  is  mainly  responsible  for  inducing  the  upward  and 
downward  ion  drifts.  When  the  wind  blows  toward  the  poles, 
a  downward  ion  drift  is  induced;  and  when  it  blows  toward 
the  equator,  an  upward  ion  drift  is  induced.  If  the  wind 
changes  direction  with  the  altitude  (a  helical  wind  shear),  the 
plasma  will  either  diverge  and  decrease  its  density  or  con¬ 
verge  and  increase  its  density  (layer  formation).  When  a  null 


in  the  wind  shear  moves  down  in  altitude,  the  ion  con¬ 
vergence  region,  and  hence,  the  intermediate  layer,  also 
descend.  Many  "active  experiments’  have  been  conducted  to 
excite  plasma  instabilities  by  creating  artificial  plasma  density 
gradients  and  observing  the  resulting  layered  structures  and 
their  nonequilibrium  dynamics.  The  active  technique  is  to 
inject  large  amounts  of  tracer  (barium  gas)  into  the  ionosphere 
from  a  rocket.  The  barium  is  ionized  by  sunlight  and,  if 
released  at  sunset  or  sunrise,  results  in  a  visible  long-lived 
plasma  made  usable  by  the  resonant  scattering  of  sunlight. 
Eigure  6.17,  from  [50],  shows  visualization  of  the  E-layer 
plasma  instabilities  using  this  Lagrangian  tracer  technique. 
The  instability  triggering  plasma  turbulence  is  the  generalized 
Rayleigh-Taylor  instability. 

Our  previous  work  has  focused  on  the  upper  troposphere 
and  lower  stratosphere  (UTLS)  shear-stratified  dynamics  of  the 
neutral  atmospheric  layers  ([51,  56,  59,  60,  62,  63,  65,  66]). 
Detailed  theoretical  investigations  and  high  resolution  nested 
calculations  were  performed  to  study  the  variability  of  the  tur¬ 
bulent  length  scales  (Tatarski,  Ozmidov,  buoyancy,  shear)  and 
the  fluxes  in  the  UTLS.  The  upper  mesosphere  and  lower 
thermosphere  region  (UMLT  internal  layer),  encompassing 
80km-150km  altitudes,  poses  similar  challenges  regarding  the 
characterization  of  the  shear  stratified  dynamics  of  the  neutral 
winds.  This  transformative  layer  marks  the  transition  to  patchy 
ionospheric  turbulence,  which  is  characterized  by  strong  scin¬ 
tillations  and  anisotropies  that  are  induced  by  the  Earth’s 
magnetic  field.  Polarized  winds  throughout  the  lower  thermo¬ 
sphere  blow  at  hurricane  speeds  of  up  to  150  m  s“^  and  produce 
enormous  wind  shears.  The  UMLT  layer  is  a  low  boundary  for 
the  E/E  layers,  so  the  information  on  the  variability  of  the  scales 
and  fluxes  is  needed  for  computational  benchmarks  and  para- 
meterizations.  Most  recently,  computer  simulation  models  of 
the  Equatorial  Spread  E  (ESE)  plumes,  or  ‘bubbles,’  have  been 
developed  ([31,  41,  58,  83,  84,  99,  100]). 

Now,  we  present  three-dimensional  numerical  results  for 
the  E  and  lower  E  region  ionosphere,  coupled  with  the  neutral 
atmosphere  dynamics.  The  inclusion  of  neutral  time-depen¬ 
dent  and  turbulent  dynamics  in  the  model  allows  us  to 
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Figure  10.  Evolution  of  density  structures. 


Figure  1 1 .  Parametric  dependence  of  the  emerging  density  structures  and  fluxes  on  the  turbulent  dynamics  of  the  neutrals. 


examine  the  charge-neutral  interactions  over  the  full  evolution 
cycle  of  an  inertial  gravity  wave  when  the  background  flow 
spins  up  from  rest,  saturates  and  eventually  breaks.  Using 
Lagrangian  analyses,  we  show  the  mixing  patterns  of  the 
ionospheric  responses  and  the  formation  of  ionospheric  lay¬ 
ers.  The  corresponding  plasma  density  in  this  flow  develops 
complex  wave  structures  and  small-scale  patches  during  the 
gravity  wave  breaking  event.  More  results  can  be  found  in 
[58,  99,  100]. 

The  heat  flux  shows  strong  vertical  mixing  in  figure  9. 
The  evolution  of  the  density  structures/interfaces  is  high¬ 
lighted  in  figure  10.  The  parametric  dependence  of  the 
emerging  density  structures  and  fluxes  on  the  turbulent 


dynamics  of  the  neutrals  is  illustrated  in  figure  1 1 .  The  plasma 
density  and  the  electrostatic  potential  are  spun  together  with 
the  neutral  flow  from  time  T  =  0.  We  show  these  two  fields  in 
figure  12.  The  E  layer  is  strongly  modified,  and  the  ion 
density  is  aligned  with  the  inertia  gravity  wave  phase. 
Figure  12  shows  the  plasma  response  to  the  wave  breaking, 
including  the  plasma  density  structure  development.  The 
electrostatic  potential  becomes  strongly  field-aligned. 

The  plasma  diffusivity  tensor,  given  in  [99],  is  controlled 
by  the  collison  frequency,  the  Boltzmann  constant  and  the 
plasma  temperature.  In  figure  12,  the  electrostatic  potential 
shows  a  strong  signature  of  field  alignment  throughout  the 
entire  simulated  column,  yet  some  inertia-gravity  wave 
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Figure  12.  Plasma  response  to  wave  breaking  and  plasma  density  structure  development.  The  electrostatic  potential  becomes  strongly  field- 
aligned. 


structure  can  still  be  observed  at  the  bottom  of  the  domain.  As 
the  flow  develops,  the  plasma  number  density  becomes 
strongly  undulated.  The  plasma  response  to  the  wave  break¬ 
ing  and  the  complexity  of  the  horizontal  structures  at  wave 
breaking  is  shown  in  figure  13.  Panel  a)  shows  the  neutral 
vertical  velocity.  Panel  b)  shows  the  ion  vertical  velocity.  It  is 
predominantly  along  the  wave  phase,  but  has  strong  three- 
dimensional  perturbations.  Note  that  the  velocity  structures 
for  the  vertical  components  of  the  neutral  and  ion  velocities 
are  quite  different.  Panels  c)  and  d)  show  the  ion  density  and 
the  electrostatic  potential,  correspondingly.  It  is  interesting  to 
note  that  the  small-scale  perturbations  along  the  wave  phase 
create  very  strong  signals  in  the  electrostatic  potential.  These 
signals  are  also  aligned  with  the  magnetic  field  line. 

To  better  reveal  the  mixing  geometry  of  the  plasma 
density,  we  use  Lagrangian  Coherent  Structures  (LCS)  from 
the  ion  velocity  to  extract  the  mixing  patterns  for  charged 
particles  ([100]).  Before  the  wave  break,  the  structures  are 
somewhat  trivial  in  that  they  simply  align  with  the  inertia 
gravity  waves  (IGW);  these  structures  highlight  the  shear 
between  the  wave  phases.  After  the  wave  break,  the  flow  is 
dominant  with  the  small-scale  structures,  which  may  lose 
their  coherence  during  shorter  times.  As  such,  we  only 
explore  LCS  during  the  wave  break  event.  We  focus  in  the 
region  where  the  large  overturning  of  the  neutral  velocity 
forces  the  ion  disturbance  to  break  up  into  patches.  The  clo¬ 
seness  in  time  reveals  the  evolution  of  the  density  patches  and 


the  LCS;  hence,  this  helps  to  explain  the  role  that  LCS  plays 
in  organizing  plasma  density  structures.  After  the  flow 
develops  coherent  structures,  we  use  the  Lagrangian  analysis 
to  characterize  its  topology  so  as  to  reveal  fine  scale  structures 
and  mixing  patterns  for  scalar  fields,  such  as  the  charge 
density.  LCS  are  made  up  of  distinguished  material  lines/ 
surfaces  that  attract/repel  local  trajectories  at  the  maximal 
rates.  In  recent  times,  the  finite-time  Lyapunov  exponents 
(FTLE)  have  conventionally  been  used  to  extract  the  LCS. 
Based  on  the  local  dynamics  of  the  trajectories,  the  FTLE 
offers  an  objective  description  of  how  nearby  trajectories 
stretch  with  the  background  flow.  As  a  result,  local  material 
lines  that  repel  nearby  trajectories  the  most  are  identified. 
Reciprocally,  attractors  are  found  to  be  local  maximizers  of 
the  FTLE  when  trajectories  are  integrated  in  backward-time. 
Note  that  large  separation  of  trajectories  could  also  be  due  to 
strong  shear. 

Eigure  14  shows  the  evolution  of  the  plasma  density  and 
the  LCS  for  a  cross-section.  The  Lagrangian  integration  time 
is  chosen  to  be  half  of  a  wave  period.  Since  the  structures  at 
this  time  are  predominantly  along  the  streamwise  direction, 
one  vertical  slice  is  sufficient  to  indicate  the  entire  flow 
structure.  In  panels  a)  and  d),  the  plasma  density  at  these  two 
times  is  shown.  It  is  apparent  that  plasma  density  striations 
break  up  into  small  lumps  and  shed  from  the  main  continuous 
wave  phase.  In  panels  b)  and  e),  the  forward-time  ETLE  field, 
together  with  the  plasma  density  iso-contour,  is  shown  at 
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Figure  13.  Plasma  response  to  wave  breaking:  complexity  of  the  horizontal  structures  at  wave  breaking. 
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Figure  14.  Evolution  of  plasma  density  and  the  Lagrangian  coherent  structures:  (a),  (d)  plasma  density,  (b),  (e)  forward  time  FTLE,  (c),  (f) 
backward-time  ETLE.  The  plasma  density  contours  show  a  correlation  between  the  LCS  and  the  plasma  density  (scalar). 


these  two  times.  The  most  dominant  feature  is  the  high¬ 
lighting  red  lines  along  the  IGW  wave  phase.  This  is  not 
surprising,  since  the  shear  in  the  polarized  IGW  will  create 
large  separation.  What  is  really  interesting  is  the  development 


of  the  vertical  highlighters  in  both  panels.  These  separating 
lines  divide  the  wave  phase  into  small  compartments,  which 
house  one  pair  of  high  and  low  density  patches.  This  pair 
evolves  together  in  the  compartment,  as  seen  from  panel  b)  to 
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panel  e).  The  vertical  highlighters  act  to  separate  different 
pairs  of  density  patches,  because  forward-time  repellers  serve 
as  boundaries  to  the  regions  where  the  scalar  patch  moves 
coherently.  In  other  words,  the  scalar  patches  on  the  left  and 
right  of  one  vertical  structure  are  going  to  separate  over  time. 
Panels  c)  and  f)  show  the  backward-time  FTLE  with  the 
density  iso-contour.  The  structures  are  predominantly  aligned 
with  the  IGW.  Note  that  the  backward-time  FTLE  correspond 
to  the  attractors;  hence,  these  newly  developed  vertical 
structures  act  to  hold  together  the  pairs  of  patches  identified  in 
the  forward-time  FTLE  compartments.  To  aid  in  the  inter¬ 
pretation  of  the  LCS,  which  acts  as  organizers  for  the  plasma 
blob  evolution,  we  focus  on  the  pair  of  blobs  shown  in  panels 
a)  and  d).  In  panel  a),  this  patch  is  still  connected  with  the  red 
wave  phase  at  the  bottom  right  of  the  panel.  The  emergence  of 
the  vertical  structures  in  panel  b)  indicates  that  on  the  left  and 
right  of  this  pair,  separated  by  the  solid  red  lines  in  panel  a) 
and  d),  motion  will  be  different  from  the  evolution  of  this 
pair.  By  examining  the  pair  of  patches  to  the  left  and  right,  we 
observe  that  the  rotation  of  the  pairs  is  different.  The  more 
complicated  vertical  structures  that  are  seen  in  panel  d) 
indicate  that  the  motion  will  be  even  more  complicated  in 
future  times.  In  addition  to  the  separating  structures,  the 
emergence  of  attracting  structures,  as  indicated  by  the  solid 
black  line  in  panel  d),  act  to  hold  the  pair  together;  so,  there  is 
not  much  difference  in  the  coherent  motion  within  the  pair. 
As  such,  we  conclude  that  during  the  overturning  event,  the 
plasma  flow  develops  separated  compartments,  in  which  pairs 
of  density  patches  evolve  together.  As  LCS  techniques  help  to 
identify  important  features  on  the  plasma  density  evolution, 
they  become  useful  in  the  extraction  of  plasma  fronts  and 
other  scintillation  producing  irregularities  that  may  affect 
operational  communication,  navigation  and  imaging  systems. 

5.  Summary  and  challenges 

Turbulent  hydrodynamic  mixing,  induced  by  the  Rayleigh- 
Taylor  (RT)  instabilities,  occurs  in  settings  as  varied  as 
exploding  stars  (supemovae),  inertial  confinement  fusion 
(ICE)  and  macroscopic  flows  in  fluid  dynamics,  such  as 
ionospheric  plasmas.  In  this  paper,  we  reviewed  physics- 
based  predictive  modeling  and  novel  multi-nesting  compu¬ 
tational  techniques  that  were  developed  in  order  to  char¬ 
acterize  strongly  inhomogeneous  non-Kolmogorov 
ionospheric  media.  Nested  numerical  simulations  of  iono¬ 
spheric  plasma  density  structures,  associated  with  the  non¬ 
linear  evolution  of  the  primary  and  secondary  Rayleigh- 
Taylor  instabilities  in  the  Equatorial  Spread  F  (ESF),  were 
presented.  For  the  limited  domain  and  nested  simulations,  the 
lateral  boundary  conditions  were  treated  via  implicit  relaxa¬ 
tion  techniques  that  were  applied  in  buffer  zones  in  which  the 
density  of  the  charged  particles  for  each  nest  was  relaxed  to 
equal  the  density  of  the  charged  particles  obtained  from  the 
parent  domain.  The  high  resolution  in  the  targeted  regions 
offered  by  the  nested  model  is  able  to  resolve  the  scintillation 
producing  ionospheric  irregularities  that  are  associated  with 
secondary  RT  instabilities,  characterized  by  sharp  gradients  of 


the  refractive  index  at  the  edges  of  the  mixed  regions  ([58]). 
In  [99]  and  [100],  our  studies  focused  on  the  charge-neutral 
interactions  and  the  statistics  that  are  associated  with  the 
stochastic  Lagrangian  motion.  We  examined  the  organizing 
mixing  patterns  for  the  plasma  flows,  which  occur  due  to 
polarized  gravity  wave  excitations  in  the  neutral  field,  using 
Lagrangian  coherent  structures  (LCS).  LCS  objectively  depict 
the  flow  topology,  and  the  extracted  scintillation  producing 
irregularities  indicate  a  generation  of  ionospheric  density 
gradients,  due  to  the  accumulation  of  plasma.  The  scintillation 
effects  that  are  induced  by  the  trapping  of  the  electromagnetic 
(EM)  waves  in  the  parabolic  cavities,  created  by  the  refractive 
index  gradients  along  the  propagation  paths,  were  analyzed  in 
[54],  [57]  and  [68]. 

The  ionosphere  plays  a  major  role  in  space  sciences  due 
to  its  important  influence  on  the  propagation  of  electro¬ 
magnetic  (EM)  waves  (radio  waves,  microwaves  and  lasers). 
High-impact  ionospheric  environments  can  significantly 
impact  communication,  navigation  and  imaging  systems, 
primarily  through  the  development  of  electron  density  irre¬ 
gularities  and  plasma  turbulence,  which  often  occurs  in  the 
vicinity  of  large  electron  density  gradients.  Associated  irre¬ 
gularities  and  inhomogeneous,  anisotropic,  non-Kolmogorov 
and  patchy  ionospheric  dynamics  can  have  a  spatial  range 
from  tens  of  kilometers  through  meter  scales.  A  wide  variety 
of  physical  processes  occur  on  these  disparate  scales,  and  this 
has  posed  a  considerable  challenge  to  the  goal  of  a  truly  self- 
consistent,  comprehensive  model-based  understanding  of 
irregularity  dynamics  and  morphology.  Recent  studies  have 
demonstrated  the  importance  of  a  full  3D,  high  resolution 
nested  approach  to  the  local  studies  of  limited-area  iono¬ 
spheric  environments  and  for  understanding  the  impact  of 
mesoscale  and  small-scale  ionospheric  processes  on  EM 
propagation  ([44,  45,  57,  58,  68-70,  99,  100]). 
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